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CHAPTER I 


Introduction 


1.1 Background 

Numerical models are critical to effective design and planning of engineered systems. 
Computers allow scientists and engineers to simulate the performance of engineered systems 
with far greater flexibility and less cost than can be achieved through physical experiments. 
Quantities that are difficult or impossible to measure experimentally can be accurately estimated 
in simulations. Simulated experiments can be repeated with variation of parameters to isolate 
cause-and-effect relationships that are important to improving designs. With accurate 
mathemati cal models, resources for physical models can be concentrated on prototypes that are 
thought to be particularly promising or used to verify certain critical modeling assumptions. 

The effectiveness of models depends on how well relevant physical phenomenon is 
represented mathematically. Many important problems in soil mechanics and particulate physics 
involve large discontinuous deformations, which are beyond the capabilities of numerical 
simulations based on continuum mechanics. Examples are soil plowing, pentrometers, pile 
driving, soil-tire interactions, hopper flows, mixing of powders, and mass movements by 
avalanche. Continuum formulations do not exist for such a wide range of behavior, particularly 
in the case of frictional materials such as sand [52]. In a particular situation, soil may deform as 
a solid, flow as a fluid, or behave as individual particles. All of these “phases” play important 
roles in the mechanical behavior, yet at present no single model exists that can account for these 
different modes of soil behavior. In the absence of such a model, a large and important class of 
soil mechanics problems lies beyond the reach of mathematical modeling. 


1.2 Towards a Continuum Mechanics Model for Granular Media 

Well-known methods of numerical approximation such as finite differences, finite 
elements or boundary elements have their origins in continuum descriptions of the media. 
Therefore, these methods are limited in principle by the assumptions of continuum mechanics. 
Historically, granular media have been treated, by Civil Engineers, as elastic or plastic solids 
because such descriptions fit expected design behavior. In fact, in many situations involving 
excavation, vehicle mobility, and material handling a fluid mechanics description for the soil 
may be more appropriate. 
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The principal difficulty in applying solid mechanics continuum theory to particulate 
media comes from the mathematical description of the kinematics that defines the movement of 
material “particles” within the continuum. The traditional use of continuum mechanics is 
limited to strains with continuous deformation. Figure 1.1 (a) shows the results from a finite 
element model of a soil plowing experiment. The continuous deformations required by the finite 
element formulation do not capture the kinematics of motion of the real soil system, which are 
illustrated in Figure 1.1 (b). Movement must obey compatibility relationships that preclude 
formation of slip planes and motions are restricted to an affine mapping in which each point in 
the initial configuration can be mapped into the deformed configuration as shown in Figure 1.2. 
The deformations thus defined represent topologically equivalent configurations because 
deformation from one configuration to another is continuous and a one-to-one mapping exists. A 
continuum can deform such that in the vicinity of any arbitrary point a second point can be found 
at a sufficiently small distance, s , in the initial configuration so that the points are separated by 
less £than in the deformed configuration. By contrast, when a discontinuity exists or a crack 
forms, two points can exist at the same initial location (e =0) 5 yet can be non-zero. 

In fact, the preceding description of a continuum is somewhat restrictive because finite 
discontinuities can exist provided relevant interface conditions are prescribed. For example, the 
Goodman joint element [21] used in finite element analyses is intended for such interface 
conditions. In fracture mechanics, discontinuities arise based on predefined fracture criteria that 
can be modeled as an evolving material topology. In these cases where compatibility breaks 
down, a mapping can exist between initial and deformed configuration by introducing conditions 
for the discontinuity. For example, while two points may initially lie at virtually the same 
location they may lie on different sides of some defined line; and, hence can be mapped to the 
new configuration uniquely. 

While the mapping between points on a discontinuity presents no special difficulties in 
principle, analysis of such points is complicated because under large deformations, points along 
the discontinuity that are initially in “contact” no longer interact mechanically after relative 
movements exceeding a grain size in magnitude. Thus normal and shear tractions carried by 
pairs of points initially, must be carried by different pairs after deformation. As a result, the 
constitutive equation that relates traction at the interface to displacement along the interface is 
more difficult to define because the traction no longer depends on the distance between the fixed 
material points A and B and the stress state is not a function of the deformation alone. 

Finally, there is the question of defining the discontinuity. In fracture mechanics, the 
formation and propagation of a fracture is controlled by a fracture criteria. No such criteria yet 
exist for shear bonding in soils. The continuum description breaks down at the point of shear 
bonding (e.g. Valanis and Peters) [52] leaving the problem ill-posed. To construct a suitable 
theory is needed. To create such a theory, a better view of the micro mechanics is needed. 

From the preceding discussion it may be concluded that the difficulties with describing 
soil, as a continuum are significant. The continuum description at best applies to the particles 
themselves, whereas in a particulate media discontinuities are the predominant feature. The 
behavior of the mass is controlled by the interactions between particles or groups of particles and 
not the character of the particles themselves. Thus, while the finite number of discontinuities in 
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a continuum can be dealt with through application of special interface relationships, the behavior 
of a particulate mass is completely dominated by interactions at interfaces. 
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Figure 1.1: Soil Plowing (a) Finite Element Analysis without provision tor discontinuity in 
deformation (b) Laboratory Experiment 
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Figure 1.2: Fracture Formation 
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The traditional approach in soil mechanics is to ignore the individual particles in the 
media and concentrate on the group or mass behavior of the soil grains as a homogenized 
continuum. Not withstanding the logical contradiction contained within the model, it has served 
well in practice because: 1) for many soil mechanics problems deformations are small, and 2) 
the deformations occur over spatial scales much larger than the largest individual particle. The 
present research deals with deformations for which neither of the two conditions is met. 

For large deformations the movement of soil particles can defy mapping by any notion 
of affine deformation. An example of such motion is the flow of sand in an hourglass (Figure 
1.3). Regardless of the initial proximity of particles A and B, they are decidedly disconnected in 
their final locations. For all but the smallest movements, the forces between particle pairs are 
independent of their initial locations. 



Figure 1.3: Sand Flow through an Hour Glass 


Further, the multiphase aspect of sand behavior is evident. In either the top or bottom chambers 
of the hourglass the sand exists essentially as a solid. The sand flowing from the top chamber 
behaves as a fluid. In the free-fall portion of the hourglass, the discrete nature of the particles 
becomes evident. Upon striking the pile in the bottom chamber the sand reverts to a fluid as it 
flows down the slope of the sand pile. Finally, the sand returns to a solid state as it comes to 
rest. 


In a fluid state, individual particle movements defy meaningful description. For fluids in 
general, traction do not depend on the details of individual particle motions but rather on time- 
averaged deformation rates. Accordingly, the traction represents a time-average of individual 
particle forces. 
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However, the fluid model is unsatisfactory for particulate media for two reasons. First, 
for small deformations, particles behave as solids whereby the stress state does depend on 
relative movements between adjacent grains. Thus, the fluid description is appropriate only 
when deformations are large. Secondly, the relationship between inter-particle force and 
movement is generally rate independent and therefore does not possess a simple viscosity law. 

In summary, it appears that particulate media poses three fundamental problem when 
viewed from the standpoint of a continuum. 

a. The motions of the individual particles cannot be mapped as affine deformations 
from the initial reference frame to a deformed configuration. 

b. Forces between material points are not related to their relative displacements if those 
displacements are greater than the nominal particle diameter. 

c. Interparticle forces do not depend on relative displacement rates. The non-viscous 
nature of flowing soil gives rise to multi-phase behavior. 


1.3 Discrete Element Method 

An alternative to the continuum description for soil and rock mechanics problems is an 
approach that models the material as a collection of individual particles that interact only at inter¬ 
particle contact points. The “particle model,” is a generic term for the class of simulation 
models where the discrete representation of a physical phenomenon involves the use of 
interacting particles. Particle models have successfully been applied to a wide class of problems 
in plasma physics, astrophysics, fluid dynamics, and molecular dynamics [6], [20], [24], [31], 
[35], and [40]. A particle model consists of a set of particles each of which has an individual 
collection of attributes (e.g., mass, particle position, velocity) and some constitutive 
relationships describing the interaction among particles. The particle attributes evolve according 
to the equations of motion. 

If a rate independent viscosity is considered by introducing a rate dependent viscosity such as 


Then 


f(t) = /u 0 => yield stress 


Cundall [1974] was one of the first to use particle-modeling techniques for evaluating 
soil and rock mechanics problems. He coined the term Distinct Element Method (DEM) and 
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developed computer programs for evaluating two-dimensional blocks of rock with complex 
shapes [15]. Cundall later developed the Universal Distinct Element Code (UDEC) for 
evaluating two- and three-dimensional interaction of a mixture of rock blocks that have different 
types of deformability [16]. UDEC provided capability for dealing with joint constitutive 
behavior, dynamic cracking, fluid flow and fluid pressure effects. In the late 1970's Cundall 
developed the computer program BALL to conduct research into the behavior of assemblies of 
discs and spheres [14]. Several authors have used DEMs to model granular assemblies’ [12], [3] 
Ng and Dobry [37] used a DEM to study small strain cyclic loading. Their simulation results 
agreed closely with trends found in laboratory tests of sands. Shukla and Sadd [45] used DEM 
to investigate how mechanical stress waves propagate in granular material and how they 
are influenced by media microstructure. Hopkins [25] used a variation of the DEM to model ice 
jamming in northern rivers and sea ice ridging in the Arctic. 

Sophisticated algorithms have been developed to describe the evolution of the particulate 
system, including the formation and breaking of inter-element contact [48]. Determination of 
forces between elements requires relationships to describe normal and shear interaction at the 
contacts. In some models, the individual particles can even “break” when stress conditions 
within the particle reach some critical level. Typically, soil particles have been modeled as two- 
dimensional circular rigid disks or three-dimensional spheres. Ting [50] has developed ellipse- 
based two-dimensional particles to represent the effects of contact flatness and particle 
angularity. Six-sided solids have also been used to model granular material [19]. 

The predominant disadvantage of DEMs is the enormous computational requirement for 
keeping track of all particle contact locations. At present, the maximum number of particles that 
can be feasibly handled in DEM computations is no more than a few tens of thousands of 
particles. Centrifuge scaling principles have been used to allow for modeling of full-scale 
geotechnical problems using a practical number of particles while maintaining stress similitude 
[48]. Of course, centrifuge scaling does not maintain geometric scaling because the ratio 
between structure to particle size in the model is much less than that of the actual media. Ba ant 
et al. [4] numerically found that the size of the particle simulation could influence the mode of 
failure. Another method to counter the computational requirement of DEM has been to couple 
DEM models with the finite element method (FEM) for analysis of rock mass [39], [2] and the 
boundary element method (BEM) for penetration tests in sand [27]. The ability to couple DEMs 
with FEMs or BEMs allows for the modeling of larger problems, while maintaining a minimal 
number of particles and reducing computer run time. However, for many problems, even this 
approach requires excessively large particle systems while adding the complexity of solving the 
coupled FEM/BEM problem. 

An initial goal of this dissertation is to increase the feasibility limit of a DEM by at least 
a factor of ten, a value that, at best, permits study of particulate systems at scales comparable to 
those of laboratory testes specimens. For example, a laboratory direct shear soil test of medium 
sand will contain over 200,000 particles. 
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1.4 Smoothed Particle (Hydrodynamic) Codes 

A class of models that admits the flexibility in kinematics description required for 
particulate modeling is the smoothed particle method (SPM). The SPM was developed by 
Monagahan [32], [33], [34], [35] originally for astronomical computation such as for galaxy 
formation and cratering. More recently the method has been used for analysis of penetration 
problems [29]. The method consists of smoothing the governing differential equation via a 
convolution integral in space thereby producing a weighted average form of the equation. By 
operating on the integral using integration by parts, spatial derivatives are moved to the 
weighting function (which is prescribed). Thus, the differential equations are converted to 
integral equations in which the dependent variable does not appear in any derivative term. The 
time evolution of the dependant variable proceeds in a manner similar to the DEM. Ostensibly, 
differentiability of the dependent variable is no longer a requirement because no derivatives are 
computed. Thus, discontinuities, separations, and change of phase can occur without 
introducing significant numerical difficulties. However, it should always be kept in mind that 
the constitutive Equations used to derive the differential Equations that serve as a starting point 
of the SPM are based on the notion of affine mapping between subsequent configurations. Thus, 
the discontinuous deformations simulated by the method may not be consistent with the physics 
on which the models are based. For example, a projectile cannot pierce an (ideal) elastic-plastic 
material because nothing in elastic-plastic theory allows such discontinuities to form in the 
continuum. The formation of the hole is strictly a numerical artifact of the spreading of the 
particles. The SPM will provide the key concepts used for the model proposed in later chapters. 


1.5 Objective of Research 

The objective of this research is to develop an alternative approach from which a 
mathematical theory for the mechanics of particulate matter can be formulated. This approach 
will be called a smooth soil particle system. Each smoothed soil particle will be modeled 
essentially as a point in the soil mass at which the state of the system is monitored. The particle 
carries with it an estimate of the state of the soil mass within its vicinity. The objective of the 
theoretical development is to derive the relationships that describe the interactions among the 
observation points. The key feature of the theory is that as the system is refined by adding more 
monitoring points the equations will begin to describe the behavior of a true particulate system. 
Importantly, the computational implementation is similar to that of traditional DEM one-to-one 
models so that computer technology developed for DEM can be incorporated into the smoothed 
soil particle system. 

The research is limited to dry, fairly uniform cohesionless soils. Specific objectives of 
the research can be listed as follows: 

a. To develop a three-dimensional discrete elements model and improve its 

computational efficiency to allow for large simulations. There are two issues driving 
this effort. First, regardless of the “smoothing” scheme ultimately developed, it is 
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likely that the general integration scheme will remain unchanged. The number of 
aggregate particles needed to model real problems in three dimensions will 
undoubtedly be large. Since, realistic problems will involve particle systems on the 
order of 100,000 to 1,000,000 particles; the need for improved computational 
efficiency is paramount. The second need for a large particle model is to develop an 
ability to model laboratory experiments on a one-to-one basis. The principal 
requirement of the averaged particle model is that it converges to the behavior of the 
one-to-one model, as the analysis mesh is refined. Therefore, for the averaged model 
to be useful, the one-to-one model must provide an accurate representation of the real 
particulate medium. To provide the link between the particle model representation 
and actual soil behavior, particle models must, at a minimum, be able to model small, 
simple experiments on soil. 

b. To demonstrate that simulation of laboratory experiments yields results comparable to 
real soils. Further, the behavior of these simulations should be controlled by well- 
defined material parameters. This objective must be met for discrete particle 
modeling to be applicable to real soils. This objective is important because the DEM 
makes simplistic assumptions as to the geometry and interaction of individual 
particles. Little work has been done to determine what properties of the DEM are 
required to realistically model soils. 

c. To establish an averaging scheme to convert properties local to the particles (e.g. 
mass, momentum) into continuum attributes (e.g. density, velocity gradients). An 
important issue to be resolved is how the coarseness of the particle system effects the 
properties of the averaged continuum. It is postulated that beyond a certain threshold 
number of particles, the actual number of particles participating in the average has 
little effect on the average. This principle, if verified, would put the averaging 
scheme on a logically sound, yet practical basis. 

d. To devise a computational procedure for modeling prototype-scale behavior using 
sparse smooth particle systems for which computations can be performed on existing 
computer systems. 


1.6 Dissertation Organization 

This Dissertation consists of seven chapters and an appendix. Chapter One covers the 
background and introduction of the dissertation, a review of particle modeling schemes, 
objective of the research, and the organization of the dissertation. The second chapter describes 
the development of the DEM model and efforts to improve the efficiency of the model that 
serves as a “numerical laboratory”. The third chapter takes up the general question of 
descriptive averaging of physical quantities of the particulate system. It also describes the 
development of continuum relationships derived based on “averaged” quantities. Chapter Four 
describes various laboratory experiments that were performed and the comparison of these test 
results with results from soil particle model simulations of the laboratory tests. Chapter Five 
describes the development of the numerical method used in the smooth particulate model. 
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Chapter Six describes the results from numerical simulations from the smoothed particle of 
laboratory soil experiments. Conclusions of the research and recommendations for further study 
are provided in Chapter Seven. 
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CHAPTER II 


The Soil Particle Model 


This chapter describes the basic mathematical formulations used to develop the soil 
particle model and the effort to optimize the model for use on high performance computing 
resources. Examples such as soil plowing and trap door simulations are shown in Figure 2.1. 
The primary objective of this study was to develop a model that could simulate laboratory soil 
experiments in which the number of computer particles is approximately the same as the 
number of actual soil particles in a laboratory experiment. 

By performing large particle simulations, large data sets of discrete particle information 
are obtained. These data are imperative to development of the smooth particle system. The 
principal requirement of the smooth particle model is that its results converge to the behavior of 
the one-to-one model, as the analysis mesh is refined. Therefore, for the averaged model to be 
useful, the one-to-one model must provide an accurate representation of the real particulate 
medium. To provide the link between the particle model representation and actual soil 
behavior, particle models must, at a minimum, be able to model simple experiments on soil. 

The objective was not to build the most physically realistic DEM model, but to build a 
model that allows for large simulations from which meaningful statistical data on the particulate 
mass can be obtained. Simplistic assumptions were made as to the geometry and interaction of 
individual particles, partly on the belief that only certain details of the inter-particle interaction 
actually cause identifiable effects in the emergent behavior of large particle groups. The simple 
particle model used in the present research contains these essential physical features. 

The developed model uses three-dimensional rigid spherical particles. Interparticle 
interactions are modeled by linear springs in the normal and tangential direction at the particle 
contacts. A local hysteretic damping law was developed to dissipate energy. To simplify the 
model, particle rotation was not allowed. It is believed that for many large deformation 
problems (e.g. plowing) with a large number of particles, individual particle rotation is not 
significant to the overall kinematics. The physical interpretation of the non-rotating particle 
assumption is that the particles possess infinite rotational inertia. Simulations by Ting, et al. 
[46] show that the effect of large rotational inertia is increased internal friction of the particulate 
media; and, they suggest that modeling real soil strength behavior requires that the rotational 
inertia of particles be increased significantly. Chapter four will show comparisons of computer 
simulations and laboratory soil tests and draws conclusions on when particle rotation is 
important. 
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Figure 2.1: Examples of DEM (a) Soil Plowing and (b) Trap Door Simulations 
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2.1 Stability of DEM 


Discrete element models generally use explicit integration schemes to determine the 
motion of the particles. This is because particle contacts are formed and lost during a time step 
leading to highly discontinuous non-linear behavior. 

Implicit integration schemes require difficult non-linear iteration, which involves the 
solution of very large systems of equations. Explicit integration schemes avoids iterations but 
their singularity comes at the expense of limitations in the computational time step. This 
section considers this limitation in detail. The critical time-step, t cr , required to ensure that the 
numerical scheme will remain stable during the simulation is approximated from a single 
degree-of-freedom system as: 



( 2 . 1 ) 


where m is the mass of the smallest particle and k is the contact spring stiffness of the 
smallest particle. In practice, the approximate critical time step is reduced by at least a factor 
of ten to ensure stability of the simulation [14]. 

Corkum and Ting [13] stated that if two particles are in contact for less than three time 
steps, energy produced from the contact is not necessarily conserved. However, if the contact 
lasts for three time steps or more, energy is conserved. Corkum and Ting derived relationships 
(Equations 2.2 and 2.6) to ensure that the contact existed for at least three time steps for 
particle-wall interaction and particle-particle interaction. For particle wall-interaction, the 
stability criteria is: 


2 



m m 


( 2 . 2 ) 


where c is a viscous damping coefficient. 

A line on a plot of n k versus n c can be defined that separates stable and unstable 
simulations: 


2-n k -n c =0 


where n k and n c are dimensionless terms defined as (see appendix A): 


and 


71 , 


kt^_ 

m 


(2.3) 


(2.4) 
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(2.5) 



ct 
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A similar relationship can be developed for particle-particle interactions. The stability criteria 
for particle-particle interactions with no global damping is: 


or 


2kt 2 

m 



( 2 . 6 ) 


2-2n k - 2n c = 0 (2.7) 

or 

n k -\-n c (2.8) 

The advantage of using these critical time step relationships is that the effect of the 
viscous damper is included in the calculation and provides a more accurate estimate of the true 
critical time step for the system. The viscous damper, as well as any other damping 
mechanism, stiffens the particulate system and reduces the critical time step. Equation (2.2) 
and (2.6) illustrate how the viscous damper stiffens the system and lowers the critical time step 
of a system. 

The issue of a critical time step is very important when trying to model real granular 
material. Sands consist of particles that have relatively small mass and large contact stiffness. 
This combination can drive the critical time step to such a small number that meaningful 
simulations cannot be performed. The dimensional analysis was originally performed 
(Appendix A) to determine if scaling laws could be used to increase the critical time step while 
maintainin g stress similitude. The results indicate that scaling laws do not provide any 
advantage to modeling the system. In fact, it was found that to maintain the stress-strain 
relationship for the material, the particle spring stiffness is dependent on the size of the particle. 
The smaller the particle is, the smaller the particle stiffness must become to maintain the stress- 
strain relationship of the material. Additionally, the critical time step was found to be 
proportional to the particle size. As smaller particle sizes are used, for a given material, the 
critical time step will decrease. This indicates that DEMs have a computational limit with 
respect to modeling small soil grains. 


2.2 Particle Force Calculations 

The interaction between two particles can be described by a set of springs, dashpots, 
and a frictional slider as shown in Figure 2.2. Normal and shear forces are computed between 
the two particles if the two particles are in contact. Two particles are in contact only if the 
distance, D, between their centers is less than or equal to the sum of their radii, i.e.: 
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( 2 . 9 ) 


D < R A + Rg 


where R A and R B are the radii for particles A and B. 

If two particles are in contact, the normal unit vector, n ,, and the shear unit vector,are 
determined. The normal unit vector n ,, points from the center of one particle A to the center 
of the particle B: 


B: - A, 


n,= 


D 


( 2 . 10 ) 


where A t and B : are the center point locations for particles A and B and the subscript i 
represents the coordinate direction. 


The shear vector t ,, is computed from the relative shearing velocity between the two 
contacting particles and is determined by: 


and 



( 2 . 11 ) 


V,. = V. -V, 


( 2 . 12 ) 


where v s is the relative shear velocity component of the two contacting particles, v" is the 
relative normal velocity component of the two contacting particles, and v' is the relative 
velocity of the two particles in contact, given by 

v '=vf-v* (2.13) 


Such that 


v. 



(2.14) 


2.2.1 Normal Force 

The normal force between the two particles, which is a repulsive force acting along the 
normal unit vector, is computed as: 


17 



r AB __ \l JlfAB 
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(2.15) 


where K AB )eJff is the normal effective stiffness constant for the two particles defined as: 


t/' AB _ 

A (/V)e# “ 


1 


1 1 ' 

]/■ A ts B 

*-N ™-N 


(2.16) 


K a n and K B are the normal stiffness constants for particles A and B, respectively, and ¥ is a 
penalty function that reduces particle penetrations. Because the stiffness of a particle is 
dependent on the particle size, K A will not equal K B when the contacting particles are of 

different sizes. A simple linear contact stiffness is used because it is believed that for large 
strain problems the nature of the force-displacement relationship at contacts is relatively 
unimportant. Rather, most deformation occurs due to interparticle movement’s [17], The 
penalty function was developed with two conditions in mind. First, 'F should equal zero when 
D-(R a + R b )> 0 . Secondly, if D - (R a + R B ) < 0 , then as the distance between particle 

centers becomes smaller, *F should increase. 


2.2.2 Shear Forces 

The shear force component of the particle-particle interaction is computed by 
integrating the relative shear velocity, v AB over time to provide the contact shear displacement, 

S AB , for the current time step, At .: 


Sf =vf At. (2.17) 

A trial shear force increment for the time step is then computed as 

(2-18) 

where K A s B ejff is the effective shear spring stiffness constant. 

As long as the two particles remain in contact, the shear force increment 
is added to the previously accumulated shear force for the contact, i.e., 
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where the indices L and L -1 refer to times t L and , and At = t L - t (L ^. 


The magnitude of the trial shear force is compared to a maximum shear force computed 
by a Coulomb-type friction law given by Equation 2.20. If the trial shear force exceeds this 
maximum shear force; the shear force is set equal to the maximum shear force: 

< 2 - 2 °) 


where <t>^ is the contact friction angle. 


2.3 Energy Damping 

The energy dissipation characteristics of the contact interaction law are critical for 
realistic simulation [14]. The major simplifying assumption that makes DEM computationally 
viable is that the particles are rigid and that interaction between particles occurs at the contact 
through discrete contact mechanisms. In real soil grains, energy is dissipated in the complex 
deformation process, which may include heat generated by plastic deformation, abrasion of the 
particle contact area, creation of sound and chipping or splitting of the particle. Particle 
interaction with boundaries may involve energy dissipation due to energy carried through the 
boundary materials. All of these mechanisms must be captured by the energy dissipation 
characteristics of the simple contact mechanisms. To improve energy dissipation 
characteristics, the introduction of a global damper was suggested by Cundall and Strack [14]. 
The global damper invokes a viscous force on a particle in proportion to its velocity, as if the 
particles were moving in a viscous fluid. Of course, for most simulations, the particles are not 
moving in a viscous fluid, leading to physically erroneous effects. 

In many problems, the static position of the particulate mass is of interest. In the static 
position, DEM represents a spring-mass (elastic) system with constant interparticle connectivity 
that, in the absence of energy dissipation will oscillate freely at its resonant frequency about the 
equilibrium position at amplitude on the order of the particle size. In real granular media the 
numerous mechanisms that exist to dissipate energy, not modeled by the DEM, ensure such 
extreme oscillations do not occur. Thus, a DEM must allow for the dissipation of energy to 
ensure that a static equilibrium position can be reached. The problem that immediately arises 
from integrating all energy dissipating mechanism into a single mechanisms is that the system 
either will be insufficiently damped or that other unwanted, physical effects will be introduced 
by the artificial energy dissipation law. The method used for energy dissipation should thus 
have the following criteria: 

1) the model should be calibrated from an “integrated” measurement of energy 
dissipation (e.g. a coefficient of restitution 

2) energy dissipation should be tied to particle interactions to avoid excessive 
dissipation during “free-fall” motions when coordination numbers (the average 
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number of contacts per particle) are small, as is the case for global viscous damping 
schemes 

3) motions should be strongly damped once a permanent contact is made to ensure 
static equilibrium can be obtained 

4) spurious rate dependence should not be introduced, as may happen when viscous 
damping is used either at the contact or as a global effect 

5) there must be a clear criteria for establishing the critical time step. 

Typically, DEM energy dissipation mechanisms are shear friction sliders, which 
dissipate energy via Coulomb friction; local viscous damping of shear and normal contact 
motions; and global viscous damping [13]. The problem with the use of a viscous damper is 
that it introduces an artificial rate effect into the problems that is not present in granular 
material. Studies of wave propagation [45] indicate that such physically inappropriate damping 
can give rise to significant errors. As noted above, in the case of the global damper the 
problem may be more acute because not only is spurious rate dependence introduced but also 
the effect is applied even in the absence of particle interactions. Additionally, the viscous 
damper stiffens the contact system, thus lowering the critical time step required to maintain 
numerical stability. 


2.3.1 Hysteretic Contact Law 

Shukla and Sadd [45] investigated three contact laws for DEM simulations; linear, 
nonlinear, and nonlinear hysteretic contact laws, and demonstrated how these laws effected the 
wave propagation process within a dry granular media. The linear and non-linear contact laws 
included viscous damping while the non-linear hysteretic contact law did not. These 
simulations where compared against experimental tests of disks made of aluminum material. 
Shukla and Sadd concluded that the non-linear hysteretic damping provided the necessary 
damping to control the inter-granular wave amplitude behavior and was the best match with the 
experimental data for both the wave speed and amplitude attenuation. 

In an attempt to improve the energy damping of the large particle system in the present 
research effort, a non-linear hysteretic damping law was developed and incorporated into the 
model. The hysteretic contact law is defined as a two-part curve consisting of a linear portion 
for initial loading and a non-linear curve for unloading and subsequent reloading. The particle 
force associated with the non-linear curve is computed as the product of the contact 
displacement, S , and a modified contact spring stiffness, K T . The function describing K r 
was arbitrarily selected as a power function that was related to the maximum contact 
displacement, S m , and the coefficient of restitution, e : 

f = k T S (2.21) 

and 
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where 


1 if = S m 

n= \2-e 

- otherwise 

_ e 

Figure 2.3 shows a typical linear loading curve and the subsequent non-linear curve for 
various levels of damping. The coefficient of restitution is defined as the height to which a 
particle will bounce after being dropped onto a rigid, flat surface divided by the initial drop 
height of the particle. The power n is related to the coefficient of restitution by considering 
that e is defined as: 



Contact Displacement 


Figure 2.3: Hysteretic Contact Law 




(2.23) 


e = 


W° 

W' 


where W 1 is the work associated with the contact during the initial loading. W 1 is defined by 
the initial linear loading force curve as: 


w' = X -k8l 


(2:24) 


W° is the work associated with the particle breaking the contact. W° is defined by the non¬ 
linear curve as: 


W° 



K 

n +1 


(2.25) 


Substituting Equations (2.24) and (2.25) into Equation (2.23) yields a relationship between the 
coefficient of restitution and n : 


2 

e =- 

w + 1 


or 


2-e 
n = - 


e 


(2.26) 


If k T at 5 m satisfies 



(2.27) 


an in view of the critical time step 


/ = 



(2.28) 


then 


t 

T c 




(2.29) 


Equation (2.29) can be used to evaluate how a change in e will affect the critical time step. 
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2.3.2 Variable Viscous Local Damper 

The hysteretic damping law will remove energy from the system only if contacts are 
formed and broken. If two particles remain in contact, they will oscillate along the non-linear 
portion of the hysteretic curve. Thus, an additional mechanism is required to ensure that the 
interacting particles will come to rest. The option of creating a more general hysteretic law is 
rejected to avoid introducing a large set of interval variables. Therefore, a variable viscous 
local damper is applied to the non-linear portion of the hysteretic damping curve. The viscous 
damper is set so that the maximum amount of viscous damping can occur without a decrease in 
the critical time step. If Equation (2.6) is used as the stability criteria, then 

m-k T t 2 -ct = 0 (2.30) 

Solving for the c that just meets this criterion: 


c = 


m-k T t 2 


Substituting Equations (2.22) and (2.27) into Equation (2.31) produces: 


1 -- 

n 


f c* "N 


c-m- 


) 


*-l 


(2.31) 


(2.32) 


or 


c = \fnkm 


1 -- 

n 


f s\ 


\ S m J 


«-l> 

) 


(2.33) 


Note that the viscous damper is “engaged” only after the contact force begins to oscillate and 
make a minor contribution to the total energy dispersion in the system. 


2.4 Particle Motion 

Once all of the contact forces acting on a particle have been determined for each 
particle, the forces are resolved into orthogonal components. A gravitational force is applied to 
each particle. The magnitude of the gravitational force is equal to the product of the mass of 
the particle and the gravitational acceleration constant for the system. Forces acting upon each 
particle are vectors summed and the instantaneous acceleration of the particle, over the current 
time step, is determined using Newton's second law of motion: 
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(2.34) 


Finally, the updated velocity 
motion: 


X? = 


2X 


m 


and location of each particle are determined using equations 


K =Jr,Vo +4ar <V o 


of 


(2.35) 


and 




i,N 


= *,\(jV-l) + At ^f(N-\) 




(2.36) 


2.5 Optimization of Particle Code 

This section describes the effort to optimize the model for use with single processor 
computing resources. Efforts to optimize the model for parallel computing resources are 
described in [8]. The preliminary modeling effort of the soil particle model focused on the 
implementation of the DEM algorithm. Little or no attention was given to computational 
performance or efficiency. The algorithm consisted of two main sections: calculation of forces 
and integration of equations of motion to move particles. Each section loops over the entire 
number of particles, performing all calculations for that section on a particle-by-particle basis. 
Contacts are determined by calculating distances between all particles, and redundant contact 
and force calculations were eliminated by identifying each particle by a unique “id” number 
and limiting the checks to the current particle and those with a larger particle “id” number. 

The determination of particle contacts consumed the bulk of the computations. The resulting 
model, though simple, could only reasonably run problems in the couple of thousand-particle 
range for any simulation requiring a large number of time steps. 

The initial performance evaluation concentrated on single processor performance on a 
CRAY Y-MP and was aimed solely at code structure. With no major algorithm modifications 
being made, the emphasis was on improved vectorization, improved input/output, and the 
elimination of redundant work. Though performance was improved by at least an order of 
magnitude, it was still insufficient. 

Performance profiling had clearly established the contact checks as the primary 
computational bottleneck. The all-against-all particle distance check, though simple, scales 
computationally as the square of the number of particles 0(N 2 ), making it impractical for large 
problems. By using a link-cell type method [23], [24] to create a neighbors list [53], the 
contact check was reduced to an O(N) operation, dramatically reducing the computational 
requirement. This method divides the physical space of the simulation into a regular grid of 
cells. 
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A particle's “owning” cell can be determined with one pass through the particles. With the cell 
size being set greater than or equal to the maximum particle diameter, possible contacts are 
limited to particles within a cell and the 26 surrounding cells. The actual implementation only 
puts particles into the upper (or lower) 13 surrounding cells to avoid duplicating redundant 
(equal but opposite) force calculations. Larger cell sizes can decrease the frequency with which 
one has to update the neighbor's list, but add considerably to the number of potential contacts. 
Since the neighbors list algorithm amounted to only 596-15% of the total runtime, it was found 
to be more efficient to minimize the cell size and update each time step. The cell size was set 
to just larger than the diameter of the largest particle in the system. This ensured that the 
number of potential contacts within a cell was minimized and that a particle's interaction was 
limited to the cell it occupied and the surrounding 26 cells. Updating each time step also 
simplifies the transition to a parallel model. 



Figure 2.4: Single processor CRAY Y-MP CPU time-per-time step vs. number of particles for 
the three phases of single processor optimization 


Figure 2.4 shows the single processor improvements for the three versions of the model 
on a CRAY Y-MP as the number of particle increases. The times given are an average time- 
per-time step. 
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Chapter III 

Continuum Representation of Discrete Particle Systems 


3.1 Background 

Most numerical models are discrete approximations to systems of differential equations. 
The differential equations are expressed in terms of continuous variables and describe the 
behavior of a continuous medium. In reality, all media are discrete and the continuum 
assumption serves only as a useful approximation that simplifies analysis. Normally, the 
difference in scale between the smallest element of the system and the problem domain is so 
great that the continuum approximation is valid. 

In the course of devising models of natural systems, two steps are essential: 1) a media 
is conceived that constitutes a continuum for which relevant balance laws can be stated and 2) 
constitutive laws are defined to close the system of equations created from the balance laws. 

The central issue of a continuum description, versus a discrete description, is the continuity of 
field quantities. In particular, descriptions of deformable continuum bodies begin with a 
requirement for deformations to be compatible, meaning that the deformation fields are 
integratible so that unique displacements can be obtained from them. Physically, compatibility 
implies that a particular set of material points can be mapped uniquely, by affine 
transformation, from one stage of deformation to the next. Such deformations do not admit 
rips, tears, and fractures or so forth within the body without additional specification of auxiliary 
discontinuity conditions. In this regard, the size of the elementary particles making up the 
medium are not an issue because even in a discrete system, individual particles can be 
constrained to move in accordance with a continuous affine mapping. For movements that 
violate affine mappings, material points must belong to independent sets, regardless of their 
size. 


In other applications, the size of the discrete elements becomes an issue because certain 
field quantities display correlation lengths whereby values of a quantity at one location are 
related to that quantity at another. Quantities associated with a discrete element display 
spatially related behavior over dimensions equal to those of the element. In media consisting of 
discrete elements, therefore, it is sometimes expedient to forego devising a continuum 
description and proceed with a model of the discrete media such as described in the previous 
chapter. Discrete models provide a basis for a more direct link to the underlying physics but 
suffer from the drawback that the numerical representation of the media is significantly coarser 
than the actual media. This limitation must be weighed against the reality of a numerical model 
of a continuum which itself is an approximation of a discrete system. It can be expected that 
any numerical approximation will be coarser than the actual media whether based on a discrete 
description or a continuum representation. 

Formulations in continuum mechanics are based on field quantities such as density, 
stress, and strain, all quantities associated with a representative elemental volume (REV). The 
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REV is taken to be that volume for which a quantity can be averaged to produce a statistically 
stationary result. A typical example is density, which is defined as the mass contained within a 
sampling volume. At the smallest scale, if the sampling volume lies within a solid grain the 
density will be that of the solid phase. If the sampling volume lies within a pore, the density 
will be zero. If an average is made for a somewhat larger volume the density will vary between 
the solid’s density and zero. As the sampling volume increases the average density will 
approach the material’s bulk density. The volume needed to obtain a statistically correct bulk 
density is the REV. For the REV formulation of the continuum, each material point is assumed 
to possess the attributes of the REV surrounding it. While it is recognized that the REV has a 
finite size, a fact that may affect accuracy of the continuum model, the REV size itself is 
typically not an explicit part of the continuum model. There is no REV size parameter. 

While it may appear that a discrete model permits a more direct link to physical 
phenomenon without incurring errors significantly greater than discrete numerical 
approximations of continuum equations, relatively few models are actually based on discrete 
models. The reason for this is simply the lack of recognition of coarse discrete models as viable 
alternatives to continuum models. Discrete models have found applications in studies of 
fundamental physical mechanisms at the fundamental scales [1], [10], [11], [28]. What is 
lacking is a general treatment of discrete models that systematically distinguishes length scales 
that arise as part of the numerical approximation process from the scales intrinsic to the 
physical system. 

Two types of discrete length scales are distinguished in the formulation that follows; a 
smoothing length scale and the characteristic physical length scale (e.g. particle size) of the 
discrete system. The numerical smoothing length scale is tied to a smoothing operator that 
distributes conserved quantities (e.g. momentum, mass, and energy) identified at a finite set of 
points within the sampling area over a mathematical continuum. The smoothing operator is 
applied to the difference equations, derived from the discrete system, that contains the natural 
length scale. For most problems, the smoothing scale is larger than the natural scale. 

However, as the smoothing scale is made smaller, the difference equations describing the actual 
discrete system are recovered. 


3.2 Smoothing Field Quantities and Strain 

There are several motivations for using a smoothing operator in lieu of a REV. First, 
the method provides a systematic way to derive equations in terms of continuous functions 
without introducing the concept of the representative elementary volume (REV) which in mm is 
tied to the properties of an infinitesimal element. In the REV approach, the characteristic 
length scales are lost. In the smoothing approach, the discrete approximation is merely 
coarsened. Secondly, the smoothed media approach provides a means to convert the statistical 
properties of the individual components, such as grain size distribution, of the discrete system 
to continuous properties. For example, properties of elastic elements in a discrete particle can 
be related to the elastic stiffness tensor. Thirdly, the smoothed continuum properties converge 
to the discrete media properties as the smoothing length goes to zero. 
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3.2.1 Implied Average Volume 

The use of a smoothing operator for modeling discrete systems has been demonstrated 
in Smoothed Particle Hydrodynamics [7], [33], [34], and [35], This averaging approach yields 
an “implied averaging volume” (IAV). The averaging is implied because the definite averaging 
volume of the REV is replaced by a space in which weighted averages are taken. For example, 
the implied mass average over the entire space is given by the integral 

m(x) = \*_j(x-x')m(x')dx' (3.1) 

where x is the location of the averaging point and x' is a location of a particle in the vicinity 
of x . The kernel function, <p(x - x'), weights the average with distance from the material 
point to which the average is assigned, giving rise to a characteristic length similar in concept 
to the REV dimension. The characteristic length, h , of the system as shown in Figure 3.1 is 
defined as the size at which the average value of a quantity becomes statistically stationary. The 
kernel function is symmetric in its argument, monotonically decays with distance, and satisfies 

f Mx-x')dx’ = \. (3.2) 

J -00 

The integral has compact support meaning 

tm ( x )= J ^(x - x')m{x') dx' (3.3) 

where 

f <f>{x - x') dx' -1 (3.4) 

J x-h 

and 

^(x-x )= 0 v{ x-x >|/i|j (3.5) 

Within the interval |x - x'| < \h\ , the kemnel (f>{x - x') is assumed to be sufficiently smooth to 
permit mass distributions of the form 
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(3.6) 


m{x)='^ j m k S(x-x') 

k= 1 

where m k is the mass of the k' h particle located at x k in the system. The Dirac delta 
operator, j(x -x* ), satisfies Equation (3.2) and is zero for x * x k . N is the total number of 
particles in the system. Equation (3.6) holds for any x including x'. Thus by substituting 
Equation (3.6) into (3.1): 

m{x) = <t>(x-x')f j m k s{x'-x k )dx'. (3.7) 

k -1 

The redistribution from discrete particle masses to the continuum of material points is mass 
conserving for the whole domain as follows 

f m(x)dx= f* r i,(x-x’)j\m k s(x'-x k )dx‘ dx. (3.8) 

J -co J -oo J x-h , , 

K=1 

Upon c hanging the order of summation and integration over x' and invoking the selector 
property of the delta operator we get 

[ m(x)dx= [ ^y\m k (j){x — x k ^dx. (3.9) 

J -oc J —oo * 

k =1 

Again, interchanging the order of summation and integration and noting the infinite extent of 
the integral, we evaluate the integral using Equation (3.2) to get 

| m{x)dx = y'^m k . (3.10) 


The effect of mass averaging is to distribute (or numerically diffuse) the particle masses 
in accordance with the weighting kernel to produce a continuous field m(x). Like the REV, 
the IAV produces a continuum in which each material point within the medium is given a 
property corresponding to the smoothed average for particles making up the (actual) discrete 
media. The key advantage of the IAV over the REV is that a characteristic size of the 
smoothing function can be incorporated into the continuum definition. 
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Figure 3.1: Characteristic Length of the Kernel Function 


3.2.2 Computation of Gradients 

A major advantage of the smoothing operator, Equation (3.1), is that derivatives and 
gradient operators for the smoothed quantities can be computed. The average gradient of the 
mass is defined as 

Vm(x)= J" + ^(x-x')Vw(x')dxr'. (3.11) 

T his formula has limited application in view of the character of m(x) given by 
Equation (3.6). However, in view of the smoothness requirements imposed on <f>(x ), 
integration by parts may be applied to the integral to obtain 

Vm(x)= J ^-V0(x-x')/n(x')dx'. (3.12) 

Note that the boundary term resulting from the integration by parts is zero when evaluated at 
the limits of x + h and x - h. Thus the gradient of the mass average becomes the weighted 
average of the mass using -V0(x - x') as the weighting kernel. The gradient is taken over 
the same volume of mass as the average itself. 



Discrete Particle System 
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3.2.3 Momentum, Velocity, and Velocity Gradient 

It is often useful to define averages in terms of weighted quantities where the weighting 
is in terms of a quantity, such as mass, that measures the attributes of the sampling space. A 
physically motivated example of a mass-weighted quantity is the average velocity, v ( (x), which 

is determined from the averaged momentum, p, (x) and average mass, m{x ): 



Pi( x ) 

m(x) 


(3.13) 


where 


p i { x )= f $(x-x , )m(x')v i (x')dx' (3.14) 

J —O0 


and m(x) is computed from Equation (3.1). 

The spatial velocity gradient, L y , is given by 

L, -W,(*)=4^yp,(*)+ z 4 T ?,(x)7m(*). (3.15) 

m{x) m\x) 

The second term in Equation (3.15) accounts for the non-uniform distribution of particles since 
m(x) varies with x . 


3.2. 4 Sampling Size Effects 

One of the primary advantages of the descriptive averaging is that size effects are 
accounted for through the use of the characteristic length of the kernel function. An 
understanding of how the characteristic length or sampling size affects the descriptive averaging 
process is required. If the characteristic length is set too small, the results from the descriptive 
averaging process will not properly model the macroscopic behavior. Conversely, if the 
characteristic length is too large, the results from the averaging process will tend to smooth out 
any localized effects, such as shear banding. 

An additional concern about the characteristic length is its effect on the numerical 
accuracy of gradients. Since gradients reflect change in a property or behavior, they require 
more data than a point variable, such as mass, to be reliable. To illustrate this point, the 
particle data set from the Treasure Island sample, described in Chapter 4, was analyzed. Each 
particle was assigned a horizontal velocity equal to the particle's vertical position. By 
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definition, this should yield a value of one for the shear component of the spatial velocity 
gradient, L n , as defined by Equation (3.15). However, some deviations from this ideal should 
be expected because in the averaging scheme the mass of a particle is lumped at the center of 
the particle. This can lead to an uneven distribution of masses that in turn will effect the 
estimate of the mass gradient. 

The shear component of the spatial velocity gradient was sampled using four different 
sampling lengths for two different kernel functions. The first kernel function is a hat function 
shown in Figure 3.1. The second kernel is a spline function [35] used in smoothed particle 
hydrodynamics. Figure 3.2 demonstrates the impact that sampling size (i.e., number of 
particles in the sample) has on the accuracy of the estimated gradient. Clearly, a fairly large 
number of particles, 200 or more is needed to deliver a good estimate of the gradient for this 
soil. 



Number of Observations 


Figure 3.2: Effects of Sampling Size on Gradient Estimate 

The spline kernel function provided a slightly better estimate of the gradient than the 
hat kernel function. However, because of the complexity of the spline kernel function and 
additional computational requirements for computing averages, the hat kernel function was used 
for the remainder of the research. 
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3.2.5 Strain 


The ability to compute the spatial velocity gradient from the discrete system allows for 
the measurement of strain in the system. Typically, Lagrangian strain is used to describe the 
str ainin g Lagrangian strain is measured relative to an initial condition. A transformation 
between the spatial coordinate grid, from where the displacement measurements are taken, and 
the material grid, Figure 3.3, which relates the system relative to an initial configuration, must 
be developed to determine the Lagrangian strain history of material points during the 
simulation. To relate points between the two coordinate systems a deformation gradient tensor 
is defined as: 



Figure 3.3: Deformed Material and Spatial Grid 
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(3.16) 


F u dXj =^dXj. 

u J dXj 

The upper case subscripts signify derivatives are taken with respect to the material coordinates. 
The lower case subscript implies the direction in the spatial coordinate system, F u maps from 

the reference grid to the spatial grid. The mapping between the spatial and reference grid must 
obey certain conditions. First, the axiom of impenetrability states that no two particles can 
occupy the same location in space at a given time. This requires that the mapping produce a 
single unique location for each particle in each grid system. Additionally, any two neighboring 
points in the reference grid must remain neighbors in the spatial grid. The mapping between 
the grids must be continuous and possess continuous derivatives with respect to space and time. 
The above conditions are met if the determinant of F is non-zero. This determinant is referred 
to as the Jacobian determinant, J. 

The deformation gradient tensor is related to the spatial velocity gradient tensor, 

l = ^L = ^L^Ll = £^L.^L (3 17) 

iJ dx: dX, dx, dtdXjdx ' 

J J J J J 

Therefore, the rate at which the deformation gradient is changing can be written as: 

Ly = FF~' (3.18) 

or 

£,=V> (3.19) 

By updating the deformation gradient tensor through time, a Lagrangian strain tensor can be 
computed as: 


Eu=\(C u S u ) (3.20) 

where S u is the Kroncker delta function and C u is the Green's deformation tensor defined as: 

C u =F u F u =F r F. (3.21) 

The rate of change for the Lagrangian strain tensor is: 

E u =~(F„F,-S u ) 3.22) 

2 at 
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or 


Eu-FnFu- (3.23) 

Finally, substituting Equation (3.19) into Equation (3.23) yields: 

^-I^FfFu ( 3 - 24 ) 

3.3 Kinematics of Averaged Quantities 

At the particle level Newton's law for linear momentum gives the kinematics 
relationship: 



(3.25) 


where C k is the total number of contacts acting on particle k . In the case of a constant 
particle mass 

Pj-mVj. (3.26) 

The smoothed average of momentum is given by the weighted average of momentum of the K s 
particles in the region: 

(327) 

Ps k=l 


where 


A=f>* (3-28) 

*=1 

and 

<f> k =</>(x i -x k ) (3.29) 

where x k is the coordinate of the particle center while the smoothed average momentum p, is 
a function of the positional coordinate x,. The average momentum rate is given as: 
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(3.30) 


Pi = 



where 


i d0 k k k 

d> = — -- = Vi v,. 

dx dt J J 


(3.31) 


It follows from Equations (3.25) and (3.27) that 


A--|i 

Ps U=l 


K‘ C K 


(3.32) 


*=1 C —1 


where 


v,= 



(3.33) 


The first term of Equation (3.32) is the convective term that arises because the averaging is 
performed over a grid (possibly fixed in space) that moves at a velocity v* relative to the 
weighted average of the particles. For simplicity of presentation, an observation frame is 
selected such that v* = 0. 


By regrouping the remaining term, the sum on the right hand side, Equation (3.32) can 
be written as 


< 3 - 34 > 

Ps /=1 

where the sum is taken over L s particle contacts within the sampling region. The superscripts 1 
and 2 denote the two particles making up the /-th contact as shown in Figure 3.4. 

Regardless of which particle is labeled 1 or 2 we have 

n) =n 2 i (3.35) 

and 

(3 36) 
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It follows that, 


and 



(3.37) 

n\ = f 2 nf 

(3.38) 


It is of interest to relate the right hand side of Equation (3.34) to a gradient of stress 
with the goal of obtaining a relationship comparable to the conservation of linear momentum for 
a continuum. To this end, the particle stress is defined as: 



(3.39) 



Figure 3.4: Labeling of Particles and their Contacts 
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where R is the radius of the particle and V is the volume of the particle. It will be shown that 
a meaningful quantity is the volume-averaged particle stress: 


n s £=1 c~l 


(3.40) 


where 


K, 




(3-41) 


k =1 

The averaged stress can also be computed from the multi-particle contact quantities as 


S, + JSV* 2 //). (3.42) 


Because </> k is the only function of position, the gradient of the averaged stress is given by (for 
uniform n s ) 


Vj5, t + V//t ! n,7 2 ). (3.43) 

/«l 

Finally, we note that in view of Equation (3.37) 

V,ff, =— £(V'*' +V/« 2 )«,/,. (3.44) 

n s i=\ 

On the other hand, we can substitute Equation (3.36) into Equation (3.34) to get 

P, =-£((«*' V:l (3-45) 

Ps /=1 

The value of can be approximated as 

<p' *<!>' +R'n)V/ (3.46) 


or as 
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4> *<t> 2 + R 2 n 2 V j <f> 2 . 


Subtracting these approximations we get: 


</>' -(f) 2 /n) + R 2 V j (f> 2 n 2 . 


( 3 . 47 ) 


(3.48) 


In view of Equation (3.35), 


<f>-(f) 2 *{r}V/ + R 2 Vj<f> 2 )n 2 (3.49) 


which gives, by substitution into Equation (3.45) 


p=— £(*'V + * 2 v/k ! /<' 

Ps /=1 


(3.50) 


By comparison of Equations (3.43) and (3.50) we have: 

PsPi = «5 V Fv 


or 


From Equations (3.26) and (3.27) 


Vjf,=—Pr 

»5 


1 K, 

— 1 ,k k k 

P = —Zj m v,. . 
Ps M 


(3.51) 


(3.52) 


Thus, 


1 


V cr =—mv i 
n r 


(3.53) 


(3.54) 


The right hand side of Equation (3.54) is dimensionally equivalent to the product of a density 
(mass per unit volume) and acceleration (v,). In fact, taking m and v, to be independent 
quantities we have 


Vjov. = pv i 


(3.55) 


39 



where 


P = 


m 


Finally, from Equation (3.55) 




where, f, is the volume averaged force given by 


1 k, 

1 XT' / k rk 


Ps M 


(3.56) 


(3.57) 


(3.58) 


Thus, the particle system when smoothed follows the same momentum balance as the 
continuum. 
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CHAPTER IV 


Large Deformation Analysis using DEM 


4.1 Introduction 

A major hindrance to developing a theory for large deformations in particulate media is 
that force and kinematics relationships are impossible to measure at the particle level in three- 
dimensional granular flow. An alternative to physical measurements is to use a simulated 
particulate system from which quantities of interest can be computed. However, any continuum 
theory resulting from such simulated “measurements” would apply to the simulated media 
rather than to actual soil. This obvious point is important in view of the fact that the particulate 
simulations grossly simplify the real interactions. However, we cannot make experimental 
measurements at the particle level, interactions are based on inference from bulk behavior 
rather than direct experimental observation. Therefore, experimental verification of the 
numerical simulation is critical. 


4.2 Previous Work on DEM Validation 

There have been a limited number of reported attempts to validate DEM simulations 
with physical experiments. Rowell [43] used experimental data from Chapuis [9] to evaluate 
two-dimensional ellipse-based DEM simulations of dense and loose biaxial tests. Chapuis 
experiments consisted of assemblies of 400 to 600 cylinders, of 10-mm length. Four particle 
diameters ranging from 19.05 to 38.76 mm were used. Rowell was able to successfully model 
the cylinder's behavior in the dense condition. However, he was unable to model the loose 
condition. Sakaguchi and Ozaki [44] compared DEM simulations with experiments of 856 
cylinders of 10 mm diameter to evaluate the formation of arches (plugging) during the flow of 
granular material. Their DEM simulations yielded flow patterns that were in good agreement 
with experimental measurements. Rong, Negi, and Jofriet [42] performed numerical 
simulations and physical experiments of the flow behavior of bulk solids in bins. The 
particulate material used in the experiments consisted of hollow acrylic cylinders, 25.4 mm 
outside diameter. Experiments were conducted in sample sizes of either 195 or 780 cylinders. 
The numerically generated particle trajectories and velocities agreed well with the observed 
experimental flow behavior. The verification experiments to date have typically been 
performed on a relatively small number (less than 10,000 particles) of idealized shaped particles 
whose motions have been restricted to two-dimensions and whose size is much larger than that 
of real soils. The boundary effects in such cases are large and the details at the particle-level 
become critical. 
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4.3 Effects of Representation Grain Size Distribution 

Much work has been done to quantify the effects of particle stiffness, rotation, and 
shape on the DEM simulation [49], and [50]. Theoretical relationships have been developed to 
assign the particle contact properties such as spring stiffness and damping to reproduce bulk 
properties of the material being modeled [45], and [47]. However, very little work has been 
done to quantify the effects of using multiple particle sizes. More important is the question, for 
accurate DEM simulations, of how many particle sizes of a real soil’s grain size distribution 
curve must be used to properly simulate the behavior of the soil. Granular soils consist of 
mineral fragments of various sizes and shapes. The ratio of largest to smallest particle diameter 
can be quite large. The representation of the grain size distribution curve can have a significant 
effect on the contact statistics and the mechanical behavior of both the actual soil and the 
simulation. To examine the effects of the representation of the grain size distribution curve on 
DEM results, a numerical study was conducted that consisted of three-dimensional simulations 
of sedimentation involving up to 93,000 particles having ratios of smallest to largest diameter 
ranging from 1:1 to 9.2:1. An actual grain size distribution curve of Treasure Island sand 
(Figure 4.1) was used as the basis for the simulations. The simulations consisted of dropping 
the computer particles into a box and allowing the particulate mass to settle due to is own 
weight. The focus of the investigation was the distribution of porosity in the simulated media, 
particularly near rigid boundaries. Evaluation of the initial porosity states of the particle 
masses included variations in porosity, distribution of contacts, and tendency of the packing to 
crystallize near hard rigid boundaries. 

The Treasure Island sand was modeled using four simulations having 1, 2, 5, and 10 
particle sizes to represent the grain size distribution curve. The weight-based grain size 
distribution curve was converted to a discrete probability distribution function represented by M 
different sized particles. The probability of a certain grain size, P(D X ), occurring in a sample 
given is given by [26] 


P(D X ) = 


M 


At ** 1 g sot2*-n 

M 


(4.1) 


where D x is the grain size diameter for a particular location on the grain size distribution curve 
(e.g. D s0 ). This equation assumes that the mass of the soil is distributed equally among the M 
different particle sizes. Clearly, as M increases, the simulation more adequately represents the 
actual soil. 

For this study, a constant sample mass of one gram was used. The number of particles 
required to represent a grain size distribution increases as more particle sizes are used, as 
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shown in Table 4.1. The increase in particles with m is due to the fact that more of the smaller 
size grains are being modeled as more particle sizes are used. Figure 4.2 plots the probability 
of particle occurrence versus grain size for each simulation of the grain size distribution. For 
well-graded soils, many more particles are needed to statistically simulate the particle 
distribution than required by uniform soils. Thus, the applicability of particle methods to 
model real soil media becomes greatly limited as the grain size distribution becomes broader. 

The representation of the grain size distribution curve can have a significant effect on 
the contact statistics and the mechanical behavior of the simulation. This fact is illustrated in 
Figure 4.3, which shows the particle packing of a simulation of only one-grain size as opposed 
to a simulation that used five different grain sizes. The simulation using only a single size 
particle size produces a very uniformly structured packing, while the simulation using five 
particle sizes produces a packing that is similar to real granular soils. 

Figure 4.4 shows how various representations of the grain size distribution curve effect 
the variation of planar porosity from a rigid boundary. A key finding is that the representation 
of the grain size distribution has a significant effect on the porosity distribution near 
boundaries. The more particle sizes used to represent a particular distribution, the less 
pronounced the boundary effect would be. When the soil is represented by only one or two 
grain sizes, the tendency is for the particles to align with the rigid wall. This alignment 
propagates inward from the wall for several particle diameters. As the number of particle sizes 
is increased, the alignment becomes limited to a distance of approximately one particle diameter 
from the wall. 

The average coordination numbers (the average number of contacts per particle) for 
each particle size in the four simulations is shown in Figure 4.5. The coordination number is 
defined as the number of contacts a particle makes with its surrounding particles. For each 
individual simulation, the coordination number is proportional to the surface area of a particle. 
Another interesting point is that the coordination number increases as the number of grain sizes 
increases, especially for the larger particles. This would indicate that by modeling more of the 
fine material tighter packing could be achieved. 

For the remaining simulations it was decided that five particle sizes would be used to 
represent the gr ain size distribution curve. This provided a realistic looking packing structure, 
while keeping the simulations to a reasonable number of particles. 


4.4 Soil Plowing Simulation 

A soil plowing experiment consisting of the horizontal translation of a vertical wall 
through a uniform Ottawa 20-30 sand was performed. The experiment was configured as a 
plane-strain test whereby motion of the wall was in one plane and the sand was confined 
between rigid glass plates. The DEM simulation was three-dimensional in that, particles were 
free to move out-of-plane to the extent permitted by the boundaries. Comparison between 
simulation and experiment focused on sand deformation, velocity of individual points within the 
mass, and the total horizontal force on the plow during advance. This particular problem was 
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chosen for study because the expected large deformations, development of shear bands and 
slope instability are characteristic of a large class of large deformation problems. 


100 



Grain Size, mm 


Figure 4.1: Treasure Island Grain Size Distribution Curve 

4.5 Plowing Experiment 

The physical model consisted of a rectangular vessel (300-mm in length x 150 mm in 
height x 8.5 mm in width) with transparent glass sides. However, only a portion of the vessel 
was used, providing test dimensions as shown in Figure 4.6. The plow was attached to the 
underside of a small, four-wheeled trolley, which traveled along two rails. The sand was 
placed in the vessel using a deep-throated funnel. After the sand was placed to the desired 
height, the left end wall of the vessel was removed, thereby allowing the sand to flow out and 
form a slope at the soil's angle of repose. As the plow advances toward this wall, sand can run 
out of the vessel. Thus, a nearly constant slope angle is maintained. The trolley is displaced at 
a constant rate of 2.5 cm/sec. 
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Video images of the plowing experiments were recorded in real time and analyzed 
using a computer vision system. An image processing and analysis program called “Particle 
Tracer” [41] semi-automatically obtained the shape of the sand surface and location of the plow 
with time, and the displacement trajectories of selected sand particles during plow advance. 

To characterize the displacement field, individual sand particles are coated with a 
fluorescent dye. The experiment is recorded on video under UV light. The bright fluorescent 
tracer particles can be segmented from the other particles in the digitized images using a simple 
thresholding operation. All pixels with a grayscale value greater than some threshold value are 
marked as foreground regions and their grayscale pixel value is set equal to 255 (white). All 
other pixels are marked as black (grayscale value = 0). The plow was coated with the same 
fluorescent dye as the tracer particles and a white background was used behind the vessel. This 
enables Tracer to also determine the soil surface profile and the location of the plow in each 
image. 



Particle Diameter, mm 


Figure 4.2: Probability of Particle Size Occurrence 


Table 4.1: Number of Particles 


Number of Particle Sizes 
(M) 

1 

2 

5 

10 


Number of Particles in Sample 

28,139 

35,156 

62,526 

93,174 
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Figure 4.3: Particle Packing using (a) One Grain Size (b) Five Grain Sizes 
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Figure 4.4: Effects of Representation of Grain Size Distribution on Porosity 
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Figure 4.5: Effects of Representation of Grain Size Distribution on Coordination Number 
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4.6 Simulation Description 

The simulation consisted of 66,544 particles, which was approximately the number of 
actual particles used in the experiment. The properties used for the simulation are shown in 
Table 4.2. The normal and shear spring stiffness was selected to keep the critical time step 
within feasible computational limits. The values used were eight orders of magnitude less than 
those estimated from elastic properties of the particles (quartz). While the low particle stiffness 
is not suitable for study of particle-scale mechanisms or for dynamic computations where wave 
propagation speeds are important, the results from the soil plowing simulation indicated that 
particle stiffness has limited effect on large flow-like deformation found in the problems under 
investigation. 


Table 4.2: Simulation Properties 


Particle Shape 

spherical 

Specific Gravity 

2.65 

Contact Stiffness: 


Normal 

1.4 kg/m 

Shear 

0.4 kg/m 

Contact Friction Angle 

15 degrees 

Particle to Wall Friction 

20 degrees 

Coefficient of Restitution 

0.04 

Plow Advance 

2.5 cm/sec 

Time Step 

2.E-5 sec 


The initial creation of the DEM sample closely followed the procedures used in the 
physical experiment. An initial particle assemblage was obtained by randomly creating 
particles in accordance with Equation (4.1) using five particle sizes. Table 4.3 shows the 
particle diameters used and their probability of occurrence. The particles were then randomly 
selected for placement on a lattice with spacing between centers large enough to minimize 
initial Interparticle forces. The particles were then “rained” into the simulated rigid-wall 
container. The sample creation simulation continues until the particles achieve static 
equilibrium. The lateral constraint is then removed from the left end of the simulated test box 
causing particles to run out, creating a natural slope. 

To reduce the size of data files and to depict discrete particle data as continuum field 
variables (e.g., density, velocity and velocity gradients, and stresses), data were mapped to a 
grid as weighted averages using the smoothing technique described in the Chapter 3. After 
smoothed averages are computed for each grid location; data visualization is accomplished 
using standard finite element post processing software. 
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Table 4.3:Discrete Ottawa 20-30 Grain Size Distribution 


Percent Passing: 

Particle Diameter 

Probability of 

10 

.63 mm 

25.90 

30 

.66 mm 

22.52 

50 

.70 mm 

18.88 

70 

.72 mm 

17.35 

90 

.75 mm 

15.35 


4.7 Results and Comparisons between Experiment and Simulation 

The locations of the tracer particles at the beginning of the experiment, at 10 mm of 
plow displacement, and at 20 mm of displacement are shown in Figure 4.7. The development 
of a soil mound ahead of the advancing plow is observed. A comparison of the trajectories of 
tracer particles numbered 1, 2, 20 and 21 (as identified in Figure 4.7) to the closest particle in 
the DEM simulation for the first 20 mm of plowing is shown in Figure 4.8. 

The simulated particles and tracer particles follow a similar trajectory. However, the 
simulated particles have a smaller net displacement than the real soil particles, particularly in 
the vertical direction. This is in part due to the softness of the simulated particles. If the 
particles had greater stiffness, they would tend to move over each other rather than compress 
between each other. Analysis of contact information at 10 mm of plow movement indicated 
that 46% of the contacts of the total simulation had penetrations greater than 1 % of the particle 
diameter and 0.8% of contacts had penetrations greater than 5% of the particle diameter. The 
highest penetrations occurred near the plow face. 

Figure 4.9 shows the horizontal and vertical particle velocity components for several 
simulated-observed particle pairs located in the plow zone over the first 20-mm of plow 
displacement. This plot demonstrates that the DEM particles are typically moving with a 
velocity magnitude that is approximately 5 mm/s slower than the observed tracer particles. 

Figures 4.10 and 4.11 reveal the shape of the plowed material at 10 mm and 20-mm 
plow displacement. Since the individual particles are too numerous to show, the DEM particle 
data was smoothed. A color spectrum is used to represent the mass densities. The color map 
ranges from yellow, which represents relatively low density, to red, which represents the areas 
of highest density. The green fringe around the perimeter results from the inclusion of the void 
space outside the slope in the averaging. Figures 4.10a and 4.11a show that the simulation 
indicates a zone of high densification just in front of the plow tip. 
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The contour lines shown on the interior of the slope in Figures 4.10 and 4.11 are 
contours of particle velocity magnitude in mm/s. A comparison of the experiment and 
simulation results reveals that the complex deformation patterns observed in the experiments in 
the vicinity of the plow are captured well by the simulation. As the plow moves into the soil, a 
zone of intense shearing forms at the tip of the plow and propagates diagonally upward to form 
a shear band. This zone of shearing moves laterally with the plow, forming a moving velocity 
discontinuity. The location of the shear band is identified by the closely spaced velocity 
contours. In the region to the left and below the 10 mm/s contour, the particles are essentially 
stationary. In the region to the right and above the 23.5 mm/s contour, sand is moving as a 
deforming plug ahead of the advancing plow. The location of the experimental contours is 
virtually identical to those obtained from the DEM simulations, indicating that the DEM model 
has predicted the location, width, and the motion of particles in the advancing velocity 
discontinuity well. 

Plow force comparisons were made between experiment of the Ottawa 20-30 sand and 
the DEM simulation. An additional plowing experiment using glass beads was also used to 
compare plow forces. The glass beads had particle sizes similar to the Ottawa 20-30 sand and 
the DEM simulation. Figure 4.12 shows the resulting plowing forces from the simulation and 
the two experiments. The Ottawa 20-30 sand had the highest plowing force, while the glass 
bead experiment provided the lowest plowing force. The low plowing force obtained during 
the glass bead plowing experiment is attributed to the glass bead's great ability to roll compared 
to the Ottawa 20-30 sand. One interesting observation is that glass beads behave quite 
differently from the Ottawa 20-30 sand. Even though the Ottawa 20-30 sand is considered 
well-rounded sand, it appears that it does not have the same particle rolling characteristics as 
the glass beads. The simulation plowing force was less than that of the Ottawa 20-30 sand. 

This was attributed to low stiffness of the particles in the simulation. 

The DEM simulation particle data was smoothed to compute the spatial velocity 
gradient tensor to examine zones of shearing and rotation within the simulation. Figure 4.12 
shows the D n term of the rate deformation tensor at 1 and 2 cm of plowing. The rate 
deformation tensor was computed from the spatial velocity gradient and is defined as: 

= (4.2) 

This figure shows that a zone of high shear is occurring just in front of the plow tip. Figure 
4.14 shows the w n term of the spin tensor. The spin tensor is also computed from the spatial 
velocity gradient and is defined as: 


It is interesting to note that even though individual particle rotations were prevented, the 
particulate system still had macroscopic rotations. 

A final observation should be made regarding the shape of the slope face. In the 
experiments the slope face was linear, while in the DEM simulation it was slightly concave 
downward. This can be attributed largely to the fact that particle rotation was suppressed in the 
simulation. This suggests that simulations involving surface slope failures such as avalanches 
or hopper flow should include particle rotation. 


4.8 Trap Door Simulations 
4.8.1 Description of Experiment 

The trap door experiment consisted of filling a narrow box of dimensions 8.5mm thick 
by 50.8 mm wide by 89.3 mm high as shown in Figure 4.15 with 15 grams of dry Ottawa 20- 
30 sand. An experiment was also performed using glass beads. A small sliding gate located at 
the bottom of the box was removed, creating a small rectangular opening at the center of the 
base. This particular problem was chosen for study because the large deformations, 
development of shear bands, and slope instability are created by gravity flow, which differs 
from the failure mechanisms in the plowing simulations. Additionally, the trap door simulation 
creates zones where high velocity gradients occur with a relatively low density of particles. 

The trap door simulation allows for the evaluation of the smoothed system ability to capture 
these sharp changes in the velocity field in sparse particle systems. 


4.8.2 Comparison of Results 

The simulation consisted of 33,270 particles, which represented 15 grams of Ottawa 
20-30 sand. The soil and interface properties in the trap door simulations were the same as in 
the plowing simulations except that the particle and wall friction angles were varied to evaluate 
their impact on the results. All simulations were run to at least one second of real-time. The 
simulation in which the interparticle friction was 25 degrees and the particle-wall friction was 
22 degrees, which had the fastest flow rate, was taken out until material had flowed through the 
box. During the Ottawa 20-30-sand experiment, flow ceased at 2.6 seconds, while in the glass 
beads experiment flow ceased at 4.3 seconds. Inspection of the tracer data revealed that the 
glass beads had more horizontal movements than the Ottawa 20-30 sand and the angle of repose 
was much flatter than that of the Ottawa 20-30 sand. In all cases the DEM simulation's flow 
rate was significantly slower than that of the actual Ottawa 20-30 sand or the glass beads. Flow 
in the DEM simulation ceased at 5.8 seconds. This was expected because by preventing 
particle rotations in the simulation the overall particulate structure was more stable and flow 
though the trap door was very difficult to maintain. Figure 4.16 shows the surface profile from 
the experiments and the DEM simulation at two seconds. Clearly the simulation did not flow at 
the rate of the Ottawa 20-30 sand or the glass beads. However, the pattern of flow for the 
simulation was similar to the Ottawa 20-30 sand. Figures 4.17 and 4.18 show the evolution of 
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the D u of the rate deformation tensor. This shows horizontal movements towards the trap 
door along the top surface. Additionally, there is a zone of high horizontal acceleration just 
above the trap door opening. 

4.8.3 Conclusions 

A comparison was made between laboratory experiments involving very large 
discontinuous deformations in sand and numerical simulations using a large-scale DEM 
computation. The magnitude of the simulation provides a unique opportunity to assess the 
validity of the DEM based on experimental results. The simulation captures the behavior of a 
particulate “continuum” while the small-scale test permits a one-to-one correspondence between 
particle gradation in the simulation and the test. The agreement between the experimental and 
simulated particle motions in the plowing experiment indicates that many details not captured 
by the simplistic particle interaction model may not be relevant in statistically large assemblies. 

These simulations represent the present limiting scale at which DEM can be used on a 
one-to-one basis. Even with improved computing hardware and algorithms, simulation capacity 
can be increased by only factors of 10 to 100. Such increases still do not translate into 
appreciably larger physical experiments. Yet it has been demonstrated that a middle ground 
does exist between computationally feasible DEM simulations and small-scale experiments that 
are applicable to engineering scale problems. 
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Figure 4.7: 


Plow Experiment at (a) 0 cm Displacement, (b) 1 cm Displacement, (c) 2 cm 
Displacement 
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Figure 4.8: Comparison of Particle Traces over 2 cm of Plowing 
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Figure 4.9: Experimental and Simulation Particle Velocities 









Figure 4.10: Comparison of Velocity Magnitudes Particle during 1 cm of Plowing (a) Simulation 
(b) Experiment 








Figure 4.11: Comparison of Velocity Magnitudes Particle during 2 cm of Plowing (a) Simulation 
(b) Experiment 
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Figure 4.13: Shear Deformation Rate at (a) 1 cm of Plowing (b) 2 cm of Plowing 
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Figure 4.14: Rotation Rate at (a) 1 cm of Plowing (b) 2 cm of Plowing 
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Figure 4.15: Trap Door Problem 
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Figure 4.16: Surface Profile at 2 Seconds for (a) Ottawa 20-30 (b) Glass Beads 
(c) Simulation 
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Figure 4.17: Horizontal Velocity Gradient with Respect to y at (a) 0.5 Seconds 
(b) 1.0 Seconds (c) 1.5 Seconds (d) 2.0 Seconds 
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Figure 4.18: Horizontal Velocity Gradient with Respect to y at (a) 2.5 Seconds 
(b) 3.0 Seconds (c) 3.5 Seconds (d) 4.0 Seconds 
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CHAPTER V 


Large Deformation Constitutive Law 


In Chapter 3, an averaging process using the concept of an implied averaging volume 
was described. This averaging technique is a methodology for obtaining continuum properties 
of a discrete system. In this chapter the descriptive averaging process is applied to the particle 
data from the soil plowing simulation described in Chapter 4. From this analysis a large 
deformation constitutive law for the smoothed continuum is developed based on a statistical 
description of the DEM particulate system. 


5.1 Governing Equations for Smoothed Particle System 

By viewing the particulate system as a continuum, the equation of motion for any 
material point within the continuum is written as: 



= -pb i +p 


5v^ 

dt 


(5.1) 


where pb, is the distributed body forces (e.g., gravity). 

In a DEM simulation, forces that result from local particle contacts control the particle 
movement. In a smoothed representation, the motion of the particles does not necessarily obey 
the equation of motion. Rather the spatial average of the contact forces satisfies the equation of 
motion. Using the method of weighted residuals to average the equation of motion over a 
region, the following equation holds: 


dy = 0 (5.2) 

where M(x- x') is a positive weighting function that weights the influence of a particle, at 
location x , as a function of its distance from the averaging point x' . 

By rearranging parts of the equation and applying integration by parts [55] to the stress 
gradient, the equation may be written as: 

dx'+ f Mcr.rijdS (5.3) 

js y j 

T his form of the equation is advantageous because the stress does not have to be continuous. 

The integration by parts produces the surface integral, which can be viewed as the traction 


f dx'= [ M(x-x') 

J* dx' . J J/? 


, dV; 
Pb,+P- 


\ r M(x-x') 


aa ‘i J. k 3*1 
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boundary condition satisfied in a weighted sense. For the remainder of this discussion it will 
be assumed that sampling is being taken away from a boundary. Applying the mean-value 
theorem to the right-hand side of Equation (5.3) yields: 


M(x - x') | 


u 5v -' 

- pb , + p — L 

H ' dt 


dV 


(5.4) 


where dV represents the volume of the region. If the region is small, the value of the integral 
can be approximated as the product of the average velocity gradient and the volume of the 
region: 


M(x - x') 


- pb, + p 


dt 


AV. 


(5.5) 


A similar procedure can be performed to the left-hand side of equation (5.3) to yield: 

dM 


dx\ 


= M(x - y') 


~ pbj + P - 


dv 

dt 


(5.6) 


Equation (5.6) implies that the equation of motion for the system can be satisfied, on the 
average, by using a sampling of the particle interactions in view of 


a u = 


=X*v, 


k n) 


(5.7) 


In the DEM simulations, the average particle stress is computed by summing the k contact 
force vectors, f t , around the particle 




(5.8) 


where / is the average distance between particle centers that are in contact, V is the volume of 
the particle, and n , is the contact normal vector. 

If the coordination number, C, is known then Equation (5.8) can be written as 

6> = (5.9) 


where is the average force vector acting on the particle. 


66 



Equation (5.9) implies that the macroscopic behavior of the discrete system could be 
approximated from a sampling of the contacts. This observation suggests that it could be 
possible to replace the highly detailed discrete particle system with a simpler model that 
represented an average particle. The average particle would maintain the same average 
properties of the total discrete system. 

The development of the smoothed particle constitutive law requires that a statistical 
description of the coordination number and the fi n j be developed. 


5.2 Analysis of DEM Plowing Simulation 

5.2.1 Material Points 

Using the plowing simulation described in Chapter 4, four material points were 
established around the plow face to monitor particle contact information during the plowing. 

The material point was defined as a point in space where descriptive averaging of the 
surrounding particle data was performed. A sampling radius of 2 mm was used. The material 
point moved with the average velocity field as computed for the surrounding particles. The 
velocity of a material point is determined using Equation (3.13). These material points provide 
a time history of each state variable (i.e. stress, strain) as related to the initial configuration of 
the material. Results from the plowing simulation were used to examine the effects of sampling 
size on calculation of state variables following the procedure of Chapter 3. Figure 5.1 shows 
the initial and final location of the four material points during the first 5-mm of plowing. Each 
point is labeled based on its initial coordinates in the x and y direction. The points were labeled 
using the notation of xxyy. Where xx is the initial location of the point in mm in the x direction 
and yy is the initial location in y. The four points are 9811, 9310, 9515, 9505. 

Point 9310 was also placed in the shear zone, but farther away from the blade. Point 
9515 was placed above the shear zone, in the region where motions were predominantly rigid- 
body-motions. Point 9505 was placed below the shear zone and was established as a point that 
had little or no movement associated with it. Each of these four sampling points exhibited 
significantly different loading characteristics, and it was expected that they would provide the 
range of information required to develop a generalized macroscopic model of the granular 
system. 


Table 5.1: Sampling Sizes for Ottawa 20-30 


Sampling Radius, h Approximate 
(mm) 

1 

2 

3 

5 


Number of Particles in Sample 

30 

250 

820 

3800 
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5.2.2 Sample Size 


Only material point 9811 was used to demonstrate the effects of sampling size. Four 
characteristic kernel lengths were used in this examination. The four sampling radii and 
approximate number of particles in the samples are shown in Table 5.1. 


Lagrangian strain was computed for the material point by methods described in Section 
3.2.4, and the Cauchy stress was computed using Equation (3.40). Figure 5.2 shows the 
change in the Jacobian determinant, which is a measure of relative volume change, as a 
function of shear strain E n , for the four sampling sizes over the first 5 mm of plowing. The 
first observation is that the shear strain decreased as the sampling sized increased. This can be 
explained by observing that, as shown in Figure 5.3, for a condition of an abrupt shearing face, 
the change in horizontal velocity with respect to the vertical axis will decrease as the vertical 
length is increased. This illustrates the point that if the sampling size is set too large, localized 
failure phenomenon will be lost. The other interesting observation from Figure 5.2 is that all 
cases except the largest sampling size show the volume contracting. The three smaller 
sampling sizes all showed the same trend of volume contraction with increased shear strain. 

The largest sampling size differed from the others in that after an initial volume contraction the 
volume began to expand. It is concluded that, at the largest sampling size, the dilation zone 
above the shear zone is influencing the results. Some of the contraction is from the 
compression of the soft particles. It is expected that the particle compression would be fairly 
constant along the face of the plow. So it is expected that the effects of sampling size in this 
region would be significantly affected by particle compression. 

Figure 5.4 plots the change in Cauchy shear stress, a u , as a function of Lagrangian 
shear strain. As expected, the maximum shear stress increased as the sampling size decreased. 
The two middle-sampling sizes have similar trends. 

It is not the intent of Figure 5.4 to suggest that there is a relationship between the 
Cauchy shear stress and the Lagrangian shear strain. As shown in Figure 5.5 the Cauchy stress 
element is measured from the spatial coordinate system, while the Lagrangian shear strain is 
measured from the material coordinate system, which deforms in the view of a spatial observer. 
The use of the Lagrangian shear strain in Figure 5.4 is simply a convenient parameter to 
measure the progress of the system. 


5.2.3 Distribution of Particle Contacts and Forces 

By obtaining the spatial velocity gradient, L tj , as described in Chapter 3, the average 

deformation field can be developed for a location within the granular media. Figure 5.6 shows 
the distortion of a unit square centered around Point 9811 during the first 5-mm of plowing. It 
can be seen that the material at Point 9811 is compressing in the horizontal direction, extending 
in the vertical direction, and shearing. 
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Sampling the contact data provides meaningful macroscopic information about the 
behavior of the granular media under extremely large deformation. Figure 5.7 shows the 
distribution of contact normal forces for all the contacts in the sampling region of Point 9811 at 
5mm of plowing. From this figure, it can be concluded that the force distribution for the 
granular material exhibits a large, random component. This should be expected because the 
granular material does not move with affine motions; but, in fact its motions are nonaffine and, 
for a region to exhibit a shear, particles must slide and move over each other. 

To provide a statistical view of the data, the data were sorted in equally spaced bins, 
based on contact angles. Thirty-degree contact angle bins were arbitrarily selected for sampling 
the particle contact data, as shown in Figure 5.8. 

Figure 5.9 plots the probability of occurrence of normal contact forces sorted by 
contact angle increment bins for plowing at 5 mm. By examining the normal contact force in 
terms of the bins, a trend from the force data begins to appear as shown in Figure 5.9. The 
area under each of these curves would represent the average force for that particular bin of 
contact. Figure 5.10 shows the force data from Figure 5.9 in percentiles. It clearly illustrates 
that the largest average normal forces are occurring in the 0 to 30-degree contact region. 

Figure 5.11 shows the average contact force for each contact angle increment bin 
during plowing. From this figure, it is seen that the average force data quantitatively follows 
the trend of the smoothed deformations shown in Figure 5.6 in that the higher average forces 
are occurring in the horizontal direction where the maximum compression is occurring. 


5.2.4 Estimate of Coordination Number 

Several researchers have developed relationships between the coordination number and 
void ratio, based on experimental data. Chang et al. [10] presented several of these 
relationships. Figure 5.12 presents the coordination number as a function of void ratio for the 
results obtained from DEM simulation. Additionally, experimental data points from other 
researchers [30], [38], and [54], presented by Chang are shown in Figure 5.12. In general, the 
data obtained from the Discrete Element Model was similar to that obtained from experimental 
data. The DEM data tended to have a lower coordination number than that of the experimental 
data. 


Analysis of the DEM has suggested that the mean stress may be a better predictor of the 
coordination number. Figure 5.13 shows the relationship of average coordination number with 
the mean stress of various material points during the plowing DEM simulation. The material 
points were selected in the regions just in front of the plow. For this simulation, the average 
coordination number can be predicted reasonably well by the mean stress of the material point. 
The data from the four sampling points shown in Figure 5.13 was fitted with both linear and 
non-linear regression. The linear regression yielded the best fit. 

C = 6.41 -0.00021 a 
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Figure 5.2: Effects of Sample on Jacobian Determinant 



Figure 5.3: Effects of Sampling Size on Shear Strain 
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Figure 5.4: Effects of Sampling Size on Shear Stress 
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Figure 5.5: Spatial and Material Coordinate Systems 
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where a < 0 for compression. It is expected that the relationship between coordination 
number and mean stress would be a function of the soil grain size distribution and the initial 
porosity. Therefore, Equation 5.10 is only valid for this specific problem and should not be 
used as a general relationship. However for the plowing problem, Equation 5.10 predicted the 
average coordination number within 3.5 percent of the measured coordination number. 

In the implementation of this relationship, certain limits are placed on the equation. 
Obviously, the coordination number cannot be negative and a lower value of 6.4 was set for the 
coordination number. There are physical considerations why the coordination number cannot 
continue to increase with linear stress and an upper limit. Under the loading conditions of the 
plow experiment, the average coordination number in the DEM simulation never exceeded 7.5, 
however experimental data from Oda [39] suggest the coordination number can go much higher 
than 7.5. 


5.2.5 Average Normal Vector 

An examination of the distribution of contact angles within the plowing simulation 
revealed that the contact angles are uniformly distributed, although the forces being transmitted 
are not. Figure 5.14 shows the distribution of contact angles at different amounts of plowing 
for the material point described in the previous section. What is interesting from this figure is 
that the distribution of contact angles remains essentially uniform even over large strains. If the 
granular material moved with affine motions, it would be expected that there would develop a 
preferred contact angle; however, this is not the case. This observation has been presented by 
other DEM simulations, in particular Ng [37] who performed a two-dimensional shearing test 
in which the ratio of the principal fabric tensor components, n u / n 22 , reached an upper limit of 
1.2. This upper limit of 1.2 in a two-dimensional case would be the equivalent of having the 
principal fabric tensor components equal to 0.45 and 0.55, these values correspond to 
distortional strain of about 5 percent if it was attributed to affine deformation. The results from 
the plowing simulation yielded a ratio of the principal fabric tensor to be less than 1.2. 
However, it should be pointed out that the plowing simulation was a three-dimensional 
simulation, using five different particle sizes, and so it was easier for the particles to arrange 
themselves in a more random fashion than the two-dimensional simulation of Ng. The ultimate 
observation of both simulations is, however, that granular material will tend to reorient its 
contacts such that the distribution of contacts remains fairly uniform with direction. If the 
contact angles were to favor a specific orientation, the granular material itself would have to 
deform like rubber. However, because the granular material consists of essentially rigid 
objects, when shearing occurs and voids are formed, material will fill the voids and the uniform 
distribution of contacts is maintained. Thus, the fabric tensor varies in accordance with 
smoothed deformation at small distortions, but reaches a limiting value at a few percent strains. 
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Figure 5.6: Smooth Deformations for Sampling Point 9881 
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Figure 5.9: Distribution of Normal Force by Sampling Bins at 5 mm Plowing 


5.3 Development of the Average Particle 

Based on the fact that force data are derived when contact forces are averaged over a fixed 
contact angle space as shown in Figure 5.9, and that the distribution of contact angles is 
uniformly distributed during the simulation, one representation of the average granular system 
could be that of a spoke system as shown in Figure 5.15. If the sampling of the DEM particles is 
sufficiently large, it follows that an average particle could be represented as having equally 
spaced contacts as represented by the spoke model. The forces at the contacts can be defined 
from the deformation that occurs as a result of the smoothed velocity field as computed by the 
spatial velocity gradient. The smoothed velocity field for the average particle can be computed 
from the estimate of the current strain rate as described by Equation (3.25). The relative 
incremental displacement for a contact direction is computed as: 

<5;=A tE u tijl (5.11) 

where is the contact vector in the initial (material) configuration and l is an average distance 
between centers of particle contact pairs. 
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Figure 5.10: Distribution of Average Normal Force by Percentile 
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Figure 5.11: Average Normal Force by Contact Angle Bins During Plowing 






(5.12) 


The normal displacement acting along the direction of the contact is determined by: 

S',' =(s r ■ n)n, 

and the incremental shear displacement is the difference between the relative and normal 
displacement: 


s; = s;-s;. ( 5 . 13 ) 

The updated normal and shear contact forces for the current time, t , are 

/;(') = //' ( '- 1) + k n S'; (5.14) 

and 

//*' ) = //('"') + K,S‘ . (5.15) 

A dry granular material does not allow tension forces, so the normal force must be 
checked for tension. If the normal force is found to be in tension, both the normal and shear 
forces for the contact are set to zero. If the normal force is found to be in compression, the 
shear force is evaluated against the friction limit: 

/ 4 < /"tan (j) (5.16) 

If the shear force is set equal to the friction limit if shear force is greater than the frictional 
limit. 

5.4 Evaluation of the Smoothed Model 

Displacement data obtained from sampling the DEM results at the four material points 
were used to derive the smoothed system to compare the predicted stress of the smoothed model 
with measured stress from the DEM simulation. 

Figure 5.16 compares the probability of contact angles of the DEM simulation and the 
data of the smoothed particle model, using 180 spokes, after 5 mm of plowing. Clearly, the 
smoothed model does not maintain the same contact angle distribution as the real granular 
material. This is because the smoothed model is only modeling the average displacements (i.e. 
modeling the granular material as an affine motion system). 
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Figure 5.12: Coordination Number Vs. Void Ratio from Chang, 1989 
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Figure 5.13: Relationship of Coordination Number and Mean Stress 
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The fact that the granular system does not follow affine motions can also be observed 
by taking a closer look at the average normal forces. Each contact force was sorted into contact 
angle increments as was shown in Figure 5.8. The normal force data was lumped into six bins 
based on contact angle. The average force was computed for each of the bins. This was done 
to get a statistical view of the data. Figures 5.17 to 5.20 plots the average of the six average 
forces and the maximum and minimum value of the average forces for the four material points. 
What is of real interest is the fact that the minimum force does not go to zero. This fact could 
also be implied from figure 5.14, which showed a nearly uniform distribution of contact angles. 
The uniform distribution of contact angles would imply that the minimum values of average 
force in any region must be greater than zero for all regions because for a contact to occur there 
must be a force. However, if the system moved with affine motion, then average forces could 
go to zero when a region lends toward tension. Figures 5.21 to 5.24 plot the ratio of the 
maximum average force to the average force and the minimum average force to the average 
force. Even though the forces in the system are increasing with plow displacement, these two 
ratios seem to reach some limit as the system goes through its deformations. This could imply 
that the contact forces are being redistributed uniformly and proportionately in angle space as 
the granular system is deformed. 

From the previous observations, it is apparent that there should be some mechanism 
that accounts for the forces and contact angles generated from the non-affine motions. One way 
of modeling the effect of the non-affine motions would be to treat this system as a convection- 
diffusion system in which the physical properties being modeled are obtained from breaking the 
model into two parts based on average motion, plus a component that is a result of the deviation 
from the average motions. This type of approach is used to model many other systems, such as 
the change in pollution concentration in groundwater flow [5]. 

A good analogy of a convection-diffusion system could be drawn by using, as an 
example, the velocity of people in a shopping mall. The interaction of the people in the mall is 
similar to the particle interaction of the discrete element code, in that people avoid coming 
closer together, since people tend to maintain a minimum separation distance. The crowd 
would flow through walkways with an organized motion; but at a closer look, the random 
movements to avoid each other’s space would be apparent. Consider if you were sitting on a 
bench in the shopping mall, acting as a spatial observation point, and you had a way of 
instantaneously sampling the velocity of each person passing you in the region. You could 
develop a spatial velocity gradient for your immediate area. If you wished to predict the 
location of a person over time, you could apply the average velocity obtained from your sample 
and integrate over time to estimate the location of the individual. Obviously, after a length of 
time, the accuracy of the estimate would degrade because no one is really travelling at the 
average velocity. The movements to avoid each person's space would create deviations from 
the average motion. Likewise, in the particle system, the use of only the smoothed velocity of 
the system will provide an increasingly poorer estimate of the current condition of the system 
with time. 
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The following presentation of the convection-diffusion system was developed from a 
presentation by Bear [5], In general the convection-diffusion system describes the change in the 
value of a parameter of interest by dividing the change into two parts. The first part is described 
by an average motion, A, and the second part by a uniformly distributed random component, B. 
For a given time step the probability that the change in the quantity is greater than A+B or less 
than A-B is zero. Within the region A-B to A+B the probability of the change of the quantity 
being x is defined as p{x)— 5/2. Because the random component is distributed uniformly 
around A. The expected value of change for a time step is A and the standard deviation would 

be (7 = B / 3. After a large number of time steps the probability of change of the value can be 
described by a normal distribution: 
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Figure 5.14: Distribution of Contact Angles During Plowing 
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Figure 5.15: Conceptual Design of Spoke Model 
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(5.17) 


where EV = TA,S 2 = TB 2 /3 ,T is the number of time steps. 

Now looking at the one-dimensional convection-diffusion equation that is described as: 


dc dc d~c 
— = V — + —7 
dt dx dx 2 


(5.18) 


where - is the change associated with the average flow of the system, D is the 


dx 


dx2 


change associated with the random component of the system. 
The solution to Equation (5.18) is 


1 


c n -jAnDt 


exp 4 Dt 


(5.19) 
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A comparison of Equation (5.19) and (5.17) shows that both solutions are of the same 
form and that 


ii 

(5.20) 

S 2 = 2 Dt 

(5.21) 


To compute the value of the quantity that results from the diffusion process, a one¬ 
dimensional discrete version of Equation (5.18) can is used: 

Ac = A tD - ~ 2C +C - (5.22) 

Ax 2 

where AC the change in the property as a result of diffusion. c° represents the smooth value 
of the property at the point of interest, c and c represent the smooth value property at points 
just in front of and just in back of the point of interest. 

The diffusion coefficient, D , is obtained from plotting the standard deviation of the 
quantity of interest with time as shown in Equation (5.21). Particle forces are carried on 
contacts that have certain a certain orientation or contact angle. Without diffusion, the change 
in contact angles would result only from the smoothed or affine motions. In order to capture 
the change in contact angle that results from the non-affme motions, the deviation of the true 
contact angle as sampled from the DEM simulation from the contact angle resulting from 
smoothed deformation was measured. 

Because the deviations in forces tend to be a function of the strain of the system, the 
time variable in Equation (5.21) must contain more than just a temporal measurement of the 
system. This situation is similar to the concept of intrinsic time has been used in endochronic 
theory to better describe the state of system [51]. Intrinsic time is defined as a monotonically 
increasing scalar of strain and time. The time intergrated magnitude of the rate of deformation 
tensor was used at the intrinsic time variable for this system. Figure 5.25 show the relationship 
between the standard deviation of the average contact angle and the time intergrated magnitude 
of the rate of deformation tensor, D t] . In general, the deviation of the contact angle increases 

with an increase of the magnitude of the rate of deformation tensor. In view of Equation 
(5.21). This would imply that the diffusion model is a good analogy for describing the random 
component of forces of the granular system. The slope of a line passing through the data points 
in Figure (5.25) represents the diffusion coefficient for the system. 
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Contact Angle 


5 mm Plowing DEM Simulation 
Smoothed System after 5mm Plowing 


Figure 5.16: Distribution of Contact Angles from Smoothed Spoke System after 
5 mm of Plowing 


By adding the diffusion to the smoothed system, the smoothed system shows more 
realistic behavior. In Figure 5.26, the ratio of the maximum average force to the average force 
and the minimum average force to the average force for the DEM simulation, plots smoothed 
system without diffusion and smoothed system with diffusion. Clearly the addition of the 
diffusion component has made the force ratio of the smoothed system behave more like that of 
the DEM simulation. Figure 5.27 compares the cumulative probability of occurrence of contact 
angle of the 12 spoke smoothed system with diffusion, a 180 spoke smoothed system without 
diffusion, and the DEM simulation. From this plot, it is observed that the smoothed system now 
maintains a contact angle distribution similar to that of the real granular system. 

Figure (5.28) to (5.31) are the results from the DEM simulation and the smoothed 
system for shear stress to shear strain for the four sampling points. The smoothed system has 
good agreement with the results of the DEM simulation. For both points 9811 and 9310, the 
smoothed system produces a yield stress that is similar to the DEM simulation. For point 9515 
which is a point at which yielding has not occurred, the smoothed system matches the increase 
of shear stress with shear strain. Even the smoothed system has good agreement at point 9505 
with the DEM simulation even though this point has been subjected to only small strain. 

Figure (5.32) to (5.39) compare the p-q plots of the DEM simulation and the smoothed 
system simulation, p is defined as: 


P = 


( 7 , + CT, 


2 


(5.23) 
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It is expected that the worst fit when comparing the DEM to the smoothed system 
would be in the p-q plots because the simulations are strain driven. Therefore, the quality of the 
model predictions is based on the ability to replicate stress. The shape of the plot for point 
9811 is predicted although there are differences in detail. For point 9310, both DEM and the 
smoothed system drift up to the correct q in about the correct range of p . For point 9515, 
the stresses are both over-predicted; although the ratio q/p is correct. Hence the good 
prediction of stress ratio. This is at a low stress range and error is probably magnified. Point 
9505 has a good prediction, especially considering the small strain. This is a very different 
stress path, yet was predicted well. In general, it is doubtful that any constitutive model could 
have done as well. The results are especially amazing considering that the parameters come 
from micro-measurements. 

Figures 5.36 to 5.39 compare the development of the stress ratio for the DEM 
simulation and the Smoothed system simulation. The stress ratio from the smoothed system 
compared very well with the measured results from the DEM simulation. 

Based on these results, it appears that the smoothed model captures most of the 
significant microscopic statistics and macroscopic mechanical behavior obtained from the DEM 
simulation. The sampling points selected for evaluating the model represented a wide variant in 
imposed strain history, but in each case the smoothed system provided the correct trend. 
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Figure 5.19: Limits on Average Normal Force for Material Point 9515 
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Figure 5.20: Average Normal Force for Material Point 9505 
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Figure 5.25: Contact Angle Diffusion 
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Figure 5.26: Comparison of Force Limits for DEM and Smoothed System, Point 9811 
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Figure 5.27: Distribution of Contact Angles from Smooth System with Diffusion 
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Figure 5.28: Comparison DEM and Smoothed System Shear Stress for Point 9310 
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Figure 5.29: Comparison DEM and Smoothed System Shear Stress for Point 9811 
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Figure 5.30: Comparison DEM and Smoothed System Shear Stress for Point 9515 
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Figure 5.31: Comparison DEM and Smoothed System Shear Stress for Point 9505 
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Figure 5.32: Comparison DEM and Smoothed System Stress Path for Point 9310 
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Figure 5.33: Comparison DEM and Smoothed System Stress Path for Point 9811 
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Figure 5.34: Comparison DEM and Smoothed System Stress Path for Point 9515 










Figure 5.35: Comparison DEM and Smoothed System Stress Path for Point 9505 
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Figure 5.36: Comparison DEM and Smoothed System Stress Ratio for Point 9310 
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Figure 5.37: Comparison DEM and Smoothed System Stress Ratio for Point 9811 
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Figure 5.38: Comparison DEM and Smoothed System Stress Ratio for Point 9515 
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Figure 5.39: Comparison DEM and Smoothed System Stress Ratio for Point 9505 
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Chapter VI 


Smoothed Particle System 


The objective of this chapter is to describe one method in which the smoothed particle 
model can be implemented for solving engineering problems. While it would be a fairly easy 
task to implement the constitutive law developed in Chapter 5 into a finite difference or finite 
element model, these numerical techniques do not allow for proper representation of the 
discontinuous soil deformations associated with large deformation problems. The new 
constitutive law was integrated into a numerical method similar to the particle-in-the-cell 
method used in smoothed particle hydrodynamic codes as described by Hockney and Eastwood 
[24], The smoothed particle system consists of two parts, particles and a fixed grid, as shown 
in Figure 6.1. The particles are considered to be material points and store information about 
the current state of the system (i.e. density, stress). A fixed grid is placed over the problem 
space. The nodes of the grid are used as collection points where the particles are sampled to 
develop smoothed velocity and displacements of the system. These smoothed velocities and 
displacements are obtained by performing a spatial interpolation of the velocities and 
displacements of the surrounding particles as described in the descriptive averaging section of 
Chapter 3. The smoothed displacement field is applied back to the particles using a similar 
spatial interpolation to obtain the spatial velocity gradient tensor for the particle. The spatial 
velocity gradient tensor is used to update the stress acting at the material point. 

Particles do not interact with the other particles directly but exchange state information 
with nodes of the fixed grid. The particles derive their motions from interpolation of the 
surrounding nodal velocities. The updated stress on the particle, which is used to update the 
forces acting on the nodes, is defined by a constitutive relationship that is derived from the 
contact information obtained from DEM simulations. 

6.1 Development of the Smoothed Particle System 

To illustrate the evolution of the smoothed particle system, a one-particle example is 
presented. In this example, only gravity is acting on the particle. For this example, the 

interpolation function will be defined as <j )', where a capital superscript indicates node number 
and a lower case superscript indicates particle number. A lower case subscript indicates 
direction. 
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In this example, a single particle is passing through a grid defined by four nodes as 
shown in Figure 6.2. A bi-linear interpolation function is used to determine the contribution of 
the particle attributes on the nodes. The four interpolation factors can be written in terms of the 
cell size, h, which defines the smoothing length scale for the problem. The four interpolation 
factors are: 


and 


* 2 


= W 1 ~{ %k +y k ) h+xk y k ] 

h 

k t k k 

_ x y 

~ h 2 

k k 

= x y 
h 2 


( 6 . 1 ) 

( 6 . 2 ) 

(6.3) 






*) k 

-xjy 

h 2 


(6.4) 


Note that the interpolation function is written so that the sum of the interpolation factors 
for a particle is one: 

^ =1 (6.5) 

/=i 


This ensures that particle quantities such as mass are conserved when transferred to the nodes. 
The derivatives of the interpolation function with respect to space are: 


d<f} X 

dx 

d£_ 

dx 

d(f> z 

dx 

d£_ 

dx 



) 


( 6 . 6 ) 

(6.7) 

( 6 . 8 ) 

(6.9) 
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dy h 2 L * 

(6.10) 

d<fr 2 x k 

dy ~ h 2 

(6.11) 

d(f> 3 x* 
dy /z 2 

(6.12) 

and 


8y h 2i J 

(6.13) 

Note: 


■y 1 dfj) 1 _ y-i 3^ _ q 

m dx ^ dy 

(6.14) 

The forces on the particle due to gravity, g, are 


o 

11 

** H 

(6.15) 

and 


fy =-m g 

(6.16) 

This force is moved to the node as follows: 


ll 

(6.17) 

The convective momentum rate is allocated to the nodes as the product of the 


momentum of the particle and the gradient of the interpolation function with respect to time 

I kill k k\ 

Pi = Pit \x y ) 

(6.18) 

where 


Pi =™ k '£t I (x k y y )v! 

7=1 

(6.19) 

+’=**' v ; + a ^' v*. 
w dx x dy y 

(6.20) 
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The mass is lumped at the nodes: 


M l =m k <t>‘ (6.21) 

1 n-2 1 

p” + 2 = P- 2 + — P i At (6-22) 

Fl M 


The velocity is thus 



M 



and the node displacement is 


d, =v” + 2A t. 


The new location of the particles is interpolated from the nodal values of d { 



(6.23) 


(6-24) 


(6.25) 


6.2 Plow Simulation 

To demonstrate the application of the smoothed particle system, a plowing experiment 
similar to the one in Chapter 4 was performed. The material properties used in Chapter 5 to 
calibrate the smoothed particle model were used in this simulation. The simulation consisted of 
placing 2200 smoothed particles into a two-dimensional box with dimensions of 100 mm in the 
horizontal direction and 22 mm in the vertical direction. The plow moved at a velocity of 25 
mm per second at a depth of 10 mm. The grid had a node spacing of 5 mm. Figure 6.3 shows 
the displacements of the smoothed particle at 5 mm, 10mm, and 20 mm of plowing. The colors 
in Figure 6.3 represent the horizontal particle velocities. The color map ranges from blue, 
which represents zero velocity to red, which represents a velocity of 25 mm/second. This 
figure can be compared to the results of the laboratory plowing experiment, Figures 4.7 and 
4.10 in Chapter 4. From this figure, similarities can be seen between the smoothed particle 
system and the granular material. In general, the horizontal and vertical motion of the 
smoothed particles is similar to that of the real granular material; however, because of the 
coarseness of the computational grid, there tends to be far more smoothing of displacements in 
the smoothed particle system and the sharp shear band zone that formed in the soil experiments 
does not form in the smoothed particle system. The same experiment was performed using a 
node spacing of 2.5 mm. Figure 6.4 shows the particle locations at 5 and 10 mm of plowing, 
when using the smaller grid. As expected, as the computational grid is refined, the smoothed 
particle system will behave more like the real granular material. Figure 6.5 shows an enlarged 
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view of both the simulations. The vectors shown in this figure represent the amount of particle 
displacement relative to the particle's original position. By examining the region just below the 
plow, it can be clearly seen that there is more smearing of the displacements at the larger grid 
size. While the smaller grid size produces more realistic results, it should be pointed out that 
the critical time step for the system becomes smaller as the size of the grid is reduced. 


(o, h) (h, h) 



(o, o) (h, o) 


Figure 6.2: One Particle Example of a Smoothed Particle System 
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Figure 6.3: Smooth Particle System Plowing 


Simulation, Grid Size = 5 mm 
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Figure 6.4: Smooth Particle System Plowing Simulation, Grid Size = 2.5 mm 


105 
















Grid Size - 5.0 mm 



Grid Size = 2.5 mm 



Horizontal Velocity mm/Sec 

WOE 

25 20 


Figure 6.5: Effects of Grid Size on Particle Displacements 
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Chapter VII 


Summary and Conclusions 


In Chapter 1, four research goals were established to develop a constitutive law for 
large deformation of granular material. The first goal was to develop a three dimensional 
discrete element model and improve computational efficiency of the discrete element model to 
allow for large simulations. The need for large particle simulations was in order to develop an 
ability to model laboratory experiments on a one-to-one basis so that the discrete element model 
could be evaluated against real soils. The second goal was to demonstrate that simulation of 
laboratory experiments by the discrete element model yields results comparable to real soils. 
Once established that the discrete element model provided a reasonable model of the real 
granular material, a third goal of establishing an averaging scheme to convert properties local to 
the particles (e.g. mass, momentum) into continuum attributes (e.g. density, velocity gradients) 
was met. From this averaging scheme a new constitutive law was developed to model large 
deformation of granular material. The fourth goal was to devise a computational procedure for 
modeling prototype-scale behavior using sparse particle systems for which computations of field 
scale problems can be performed on existing computer systems. The computational procedure 
used a particle scheme to ensure proper representation of the discontinuous soil deformations 
associated with large deformation problems. 


7.1 Significant Findings 

a) DEM comparison with Soil Experiments - A comparison was made between 
laboratory experiments involving very large discontinuous deformations in sand and numerical 
simulations, using a large-scale DEM computation. The magnitude of the simulation provides a 
unique opportunity to assess the validity of the DEM, based on experimental results. The 
simulation captures the behavior of a particulate “continuum” while the small-scale test permits 
a one-to-one correspondence between particle gradation in the simulation and the test. The 
agreement between the experimental and simulated particle motions in the plowing experiment 
indicates that many fine-grained details not captured by the simplistic particle interaction model 
may not be relevant in statistically large assemblies. 

b) Continuum representation of a discrete system - An averaging process using the 
concept of an implied averaging volume was developed for obtaining continuum properties of a 
discrete system. Effects of sampling size were quantified, which lead to the conclusion that a 
significant sampling size (number of particles in the sample) is required to obtain meaningful 
macroscopic properties of the granular media. 

c) Micro-Mechanical Modeling of Granular Media - The analysis began with 
consideration of a smoothing of the DEM quantities which amounts to application of a weighted 
residual approximation of the difference equations governing the DEM simulations. The 
smoothing process eliminates spatial detail and only the statistical descriptions of particle 
interactions are needed to evaluate the equations of motion for the smoothed system. The key 
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to a continuum theory for granular media, therefore, is to relate the statistics of particle 
interactions to the kinematics of the smoothed system. This analysis implies that the 
macroscopic behavior of the discrete system could be approximated from a sampling of the 
contacts. This observation suggests that it could be possible to replace the highly detailed 
discrete particle system with a simpler model that represented an average particle. The average 
particle would maintain the same average properties of the total discrete system. It is 
concluded that without a micro-mechanical approach based on physical measurements, a 
satisfactory theory would be difficult to develop. 

d) Non-Affine Motions - The evolution of contact properties is not readily determined 
from averaged particle movements because of non-affine components of particle interaction. A 
diffusion component was added to the smoothed model to account for the non-affine motions 
that are present in the granular media. The addition of the diffusion component made the 
statistics of the smoothed system behaves more like that of the DEM simulation. 


7.2 Need for Future Research 

The research presented in this dissertation was limited to dry, uniform, granular 
materials. Additional research should be conducted to quantify the effects that particle shape 
and grain size distribution has on the smoothed particle system. The smoothed particle system 
should be evaluated against geotechnical problems different from the plowing problem, to 
verify its application to more general problems. As computational capabilities improve, 
research should be conducted to model fine grained, clay particles. 

Additionally, work should be conducted to include pore fluid into the system to allow 
for more realistic modeling of in-situ soil conditions. One of the predominant disadvantages 
with the traditional discrete element method is the inability to easily model a pore fluid. The 
inclusion of the pore fluid into the smoothed particle system should be comparable to a coupled 
deformation-fluid flow model derived from Biot's consolidation theory, because the smoothed 
particle system is defined as a continuum. 
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APPENDIX A 

Dimensional Analysis of Particle System 


A.1 Development of #-Numbers 

A n -number is a dimensionless number that is developed from a collection of key 
parameters that describe a physical process of interest. The n -number is used to identify 
proper scaling relationships for maintaining similitude between a prototype model and the full- 
scale system. This section describes a method for determining n -numbers for a system and 
has been adapted from the presentation of Freitag [18]. 

The system being evaluated can be described by a set of physical parameters, 

Pj P 2 ,• • • P n , the fundamental dimensions are d x , d 2 , • • ■ d r so that, dimensionally: 

p s S=d e x u d e 2 u ■■■d e ; s (A.l) 


where the symbol = means “dimensional equal to” and the e rs are the components indicating 
the multiplicity of each dimension in the parameter. 

For example, the parameter mass, m, is dimensionally equal to 

m=Fr l T 2 (A. 2) 

where F, L, and T represent the dimensions force, length, and time. 

The 7i -number is developed by combining the parameters in such a way as to yield a 
dimensionless number. 


A1 


jt = P?P?P? (A.3) 

To obtain a dimensional relationship, the sum of the exponents of each dimension must 
be zero. Thus, 
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(A.4) 
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* m *« 


where the number of dimensions, r, is less than the number of parameters, n. 

The values of the e tJ are known from the nature of the parameters, but the x ; must be 
found. The set of equations representing the relations among the unknowns x, is governed by 
rules of matrix algebra. The coefficients of the x,. from a matrix of particular rank, with the 
columns representing the constitution of each parameter and the rows, the frequency of each 
dimension as indicated. 
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(A. 5) 


There are n-r more unknowns than equations, so the solution consists of expressions for r of the 
unknowns, x ( , in terms of the other n-r unknowns. To ensure that the matrix solution used 
here will provide the maximum number of linearly independent solutions, the equations and the 
unknowns must be arranged to obtain a non-zero determinant in the upper left corner of the 
matrix. Standard matrix manipulation is used to obtain an r x r identity submatrix, as shown 
below: 
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(A. 6) 


A2 



Each row of the solution matrix yields an expression of one of the r parameters in terms of the 
remaining (n-r) parameters. The solution can been obtained by means of elementary 
simultaneous equations to yield the following expressions. 


*1 
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C Xn X n 
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C 2n x n 
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^~'r,r+2 X r+2 

+ 

C 2 n X n 


Assigning the values x r+l =1 and all others, x r+2 , ... x n = o for the first solution, it 
is found from Equation A.7 that x ] =• C lrM ,x 2 -C 2 r+2 ,...x r =C, r+1 . Similarly, if x r+l =0 
and x r+2 =1, all x r+2 to x n =0, it is found that x x - C lr+2 x 2 = C 2r+2 ...x r = C r r+2 . This 

process can be continued until all (n-r)x terms from x r to x„ have successively been assigned 
the value of 1. These results can be arranged in matrix form. 
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(A. 8) 


The submatrix on the right is an (n-r) identity matrix. This submatrix has a nonzero 
determinant, so the rank of the matrix is n-r, which is equal to the number of rows. 
Consequently, the rows in the matrix are linearly independent, and the rows constitute a 
fundamental system of solutions. The numerical values in each row are the set of exponents of 
the corresponding parameters that make up the n term, thus 


r r r 

— „ M.r+m „ ^3 .r+m „ 

n m=P\ Pi Pi ■■■Pr 


(A.9) 


where 1 <m<n — r. 

By inspection, the equation of n terms (Equation A.9) can be written from the matrix 
of solutions (Equation A.7). The matrix of the coefficients in their first r columns of Equation 
A.9 is simply the transpose of the r x (n-r) submatrix on the right of Equation A.7. Note, 
however, that the signs of all coefficients must change at the same time to be equivalent to the 
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operation represented by Equation A.8. This matrix is then simply augmented by the (n-r) x 
(n-r) identity to yield the form shown as Equation A.9. 


Application of n -Theorem to the Soil Particle System 

Eight parameters were chosen to describe the soil particle system. These parameters 
included stress, gravitational constant, time, density, mass, particle stiffness, particle viscous 
damping coefficient, and particle diameter. 

These 8 parameters can be recorded with their dimensions in matrix form as follows. 



G 

g 

T 

r 

m 

k 

C 

d 

F 

i 

0 

0 

i 

i 

1 

i 

0 

L 

-2 

i 

0 

-3 

-i 

-1 

-l 

1 

T 

0 

-2 

1 

0 

2 

0 

i 

0 


(A. 10) 


The determinant of the matrix is non-zero, so the rank of the matrix is three. 
Therefore, eight parameters minus the rank of three yield five n terms that will describe the 
system. By using the method described in the previous section, the five n terms can be 


derived as: 
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Use of n -numbers for Scaling Relationships 

Because a n number is dimensionless, it can be combined with other n numbers to 
produce other dimensionless relationships among the parameters. For example, dividing n, by 

tt 2 yields: 
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which by rearranging 
of the system. 
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the terms will produce the form of the equation for the critical time step 



(A. 14) 


Another interesting relationship occurs by dividing n 5 by n z to get 



(A. 15) 


This equation can be combined with the parameter strain, £ , which is a dimensionless number, 
to develop a stress-strain dimensional relationship. 


a = 


7ik 

- £ 
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(A. 16) 


where — could be defined as a Young's modulus for the media. This indicates that for a 
d 

given material, i.e. quartz, the contact spring constant must vary with the radius of the particle. 
Additionally, because the mass of a particle is proportional to the cube of the particle radius, 
equation A. 14 implies that the critical time step will become smaller as smaller particles are 
used. 


A typical use of n -numbers in geotechnical engineering is in the development of 
centrifuge scaling laws. In a centrifuge test, the prototype model is subject to a large 
gravitational force. This change, a , is represented by 


a — 


9 profype 

9 fullscale 


(A. 17) 
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Table A.l: Scaling Relationships for Centrifuge Test 

Quantity _ 

Length, d 
Contact stiffness, k 

Viscous damping,c 

Mass of particle 
Unit weight 
Mass density 
Time 
Stress 
Strain 

The objective of the centrifuge scaling laws is to maintain stress similitude between the 
prototype and the full-scale system (Le.,cr p = <7 fs ). The Table A.l summarizes the scaling 
relationship for the particle system to a centrifugal test. 

The dimensional analysis was originally performed to determine if scaling laws could 
be used to increase the critical time step while maintaining stress similitude. However, the 
results of the dimensional analysis indicate that scaling laws do not provide any advantage to 
modeling the system. 


Full Scale 


Centrifugal Model a a g's 


1 


1 la 
\la 

\la 2 

1/a 3 

a 

1 

1/a 

1 

1 


A6 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMBNo. 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 
data needed, and completing and reviewing this collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing 
this burden to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202- 
4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently 
valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. ____ 

1. REPORT"DATE (DD-MM-YYYY) 2. REPORT TYPE ~ 3. DATES COVERED (From - To) 

August 2000 _Final report___ 


4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER 

Application of DEM to Micro-Mechanical Theory for Large Deformations of Granular 5b. GRANT NUMBER 

Media _ 

5c. PROGRAM ELEMENT NUMBER 


6. AUTHOR(S) 

David A. Horner and John F. Peters 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Engineer Research and Development Center 
Geotechnical Laboratory 
3909 Halls Ferry Road 
Vicksburg, MS 39180-6199 


9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

U.S. Army Engineer Research and Development Center 
3909 Halls Ferry Road 
Vicksburg, MS 39180-6199 


5d. PROJECT NUMBER 


5e. TASK NUMBER 


5f. WORK UNIT NUMBER 


8. PERFORMING ORGANIZATION REPORT 
NUMBER 

ERDC/GL TR-00-7 


10. SPONSOR/MONITOR’S ACRONYM(S) 


11. SPONSOR/MONITOR’S REPORT 
NUMBER(S) 


12. DISTRIBUTION / AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited. 



14. ABSTRACT 

A constitutive theory is developed for granular material undergoing arbitrarily large deformations. A three-dimensional discrete 
element model (DEM) was developed to simulate granular material. The computational efficiency of the DEM was improved to allow for 
modeling of large particle systems. The need for large particle simulations was to develop on ability to model laboratory experiments on a 
one-to-one basis so that the discrete elements model could be evaluated against real soils. A comparison was made between laboratory 
experiments involving very large discontinuous deformations in sand and numerical simulations using large-scale DEM computation. The 
magnitude of the simulation provided a unique opportunity to assess the validity of the DEM, based on experimental results. The agreement 
between the experimental and simulated particle motions in the plowing experiment indicates that details not captured by the simplistic 
particle interaction model may not be relevant in statistically large assemblies. Once it was established that the discrete element method 
provided a reasonable model for real granular material, an averaging scheme was developed to convert properties local to the particles (e.g., 
mass, momentum) into continuum attributes (e.g., density, velocity gradients). From this averaging scheme, a new constitutive law was 
developed to model large deformation of granular material. 


15. SUBJECT TERMS 

Discrete element modeling. Large deformation, Numerical modeling 


16. SECURITY CLASSIFICATION OF: 


a. REPORT 

b. ABSTRACT 

C. THIS PAGE 

UNCLASSIFIED 

UNCLASSIFIED 

UNCLASSIFIED 


17. LIMITATION 
OF ABSTRACT 



19a. NAME OF RESPONSIBLE PERSON 


19b. TELEPHONE NUMBER (include area 
code) 


Standard Form 298 (Rev. 8-98) 

Prescribed by ANSI Std. 239.18 

























